Quantum criticality, lines of fixed points, and phase separation in doped 
two-dimensional quantum dimer models 
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We study phase diagrams of a class of doped quantum dimer models on the square lattice with 
ground-state wave functions whose amplitudes have the form of the Gibbs weights of a classical 
doped dimer model. In this dimer model, parallel neighboring dimers have attractive interactions, 
whereas neighboring holes either do not interact or have a repulsive interaction. We investigate the 
behavior of this system via analytic methods and by Monte Carlo simulations. At zero doping, we 
confirm the existence of a Kosterlitz-Thouless transition from a quantum critical phase to a columnar 
phase. At low hole densities we find a dimer-hole liquid phase and a columnar phase, separated by 
a phase boundary which is a line of critical points with varying exponents. We demonstrate that 
this line ends at a multicritical point where the transition becomes first order and the system phase 
separates. The first-order transition coexistence curve is shown to become unstable with respect to 
more complex inhomogeneous phases in the presence of direct hole-hole interactions. We also use a 
variational approach to determine the spectrum of low-lying density fluctuations in the dimer-hole 
fluid phase. 



I. INTRODUCTION 

The behavior of doped Mott insulators is a long- 
standing open and challenging problem in condensed- 
matter physics. Mott insulators are the parent states of 
all strongly correlated electronic systems and as such play 
a crucial role in our understanding of high-T c supercon- 
ductors (HTSC) and many other systems. Their strongly 
correlated nature implies that their behavior cannot be 
understood in terms of weakly coupled models. Except 
for the very special case of one spatial dimension, the 
physics of doped Mott insulators is currently only under- 
stood at a qualitative level. The solution of this challeng- 
ing problem remains one of the most important directions 
of research in condensed-matter physics. 

Quantum dimer models (QDM)A provide a simplified, 
and rather crude, description of the physics of a Mott in- 
sulator. They provide a correct description of the physics 
of Mott insulators in regimes in which the spin excitations 
have a large spin gap. QDMs were proposed originally 
within the context of the resonating-valence bond (RVB) 
mechanism of HTSCi^ These systems are of great in- 
terest as they can yield hints on the behavior of more 
realistic models of quantum frustration. 

The main idea behind the formulation of QDMs is that, 
if the spin gap is large, the spin degrees of freedom be- 
come confined in tightly bound singlet states which, in 
the extreme limit of a very large spin gap, extend only 
over distance scales of the order of nearest-neighbor sites 
of the lattice. Thus, in this extreme regime, the Hilbcrt 
space can be approximately identified (up to some impor- 
tant caveats^) with the coverings of the lattice by valence 
bonds or dimers. 

Surprisingly, even at the level of the oversimplified pic- 
ture offered by QDMs, the physics of doped Mott in- 



sulators remains poorly understood. In this paper we 
explore the phase diagrams, and the properties of their 
phases, of QDMs generalized to include interactions be- 
tween dimers (or valence bonds), and between dimers 
and doped (charged) holes. Basic aspects of the physics 
of these models are reviewed in Ref. and references 
therein. 

Undoped QDMs have been studied more ex- 
tensively and by now they are relatively well 
understood 1 1 1 1 On bipartite lattices, their 

ground states show either long-range crystalline valence 
bond order of different sorts or are quantum criticalf^^ 
while on non-bipartite lattices their ground states are 
typically disordered and are topological fluids^ 

The more physically relevant doped QDMs, with 
a finite density of charge carriers (holes), are 
much less understood, although some properties are 
known, 12 ' 13 ! 14 ! 15 ! 16 ' 17 ! 18 In QDMs a spin-± hole fraction- 
alizes into a bosonic holon, an excitation that carries 
charge but no spin, and a spinon, a fcrmionic excitation 
that carries spin but no charge J^ddl Holons can be re- 
garded as sites that do not belong to any dimer, whereas 
spinon pairs are broken dimers. This form of electron 
fractionalization is observable in the spectrum of these 
systems only in the topological disordered (spin-liquid) 
phases of the undoped QDM. Otherwise, as in the case 
of the valence bond crystalline states which exhibit long 
range dimer order, spinons and holons are confined and 
do not exist as independent excitations^ 

In this paper, we consider several interacting QDMs on 
a square lattice at finite hole doping, and discuss their 
possible phases and phase transitions as a function of hole 
density and strength of the interactions. At any finite 
amount of doping the system will have a finite density of 
holes, which are hard core charged bosons in this descrip- 
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tion. To simplify the problem, in this work we do not con- 
sider the physical effects of the charge-neutral fermionic 
spinous which in principle should also be present. Thus, 
at this level of approximation, all spin-carrying excita- 
tions are effectively projected out. The remaining de- 
grees of freedom are thus dimers ("spin-singlet valence 
bonds") and charged hard-core bosonic holes. Already 
this simplified picture of a strongly correlated system is 
very non-trivial. 

For a certain relation between its coupling constants, 
known as the Rokhsar-Kivelson (RK) condition,— QDM 
Hamiltonians, both with and without holes, can be writ- 
ten as a sum of projection operators. These RK Hamil- 
tonians arc manifestly positive definite operators. For 
this choice of couplings, the ground-state wave function 
is a zero-energy state which is known exactly. This RK 
wave function is a local function of the degrees of free- 
dom of the QDM, the local dimer and hole densities. The 
quantum-mechanical amplitudes of the RK wave func- 
tions turn out to have the same form as the Gibbs weights 
of a two-dimensional (2D) classical dimer problem with a 
finite density of holes. For the generalized doped QDMs 
that we consider, the norm of the exact ground-state 
wave function is equal to the partition function of a sys- 
tem of interacting classical dimers at finite hole density. 
This mapping to a 2D classical statistical mechanical sys- 
tem, for which there is a wealth of available results and 
methods, makes this class of problems solvable^^ i 11 ' 17 

In this work wc will investigate the behavior of doped 
QDMs which satisfy the RK conditions by studying the 
correlations encoded in their ground state wave functions. 
The phase diagrams of these systems turn out to be quite 
rich. As we shall see, these simple models can describe 
many aspects of the physics of interest in strongly cor- 
related systems, including a dimer-hole liquid phase (a 
Bose-Einstein condensate of holes), valence-bond crys- 
talline states, phase separation, and more general inho- 
mogencous phases. The undoped version of this system 
was studied in detail in Rcf. 1 141 , where a quantum phase 
transition was found that was argued to belong to the 
Kosterlitz-Thouless (KT) universality class, from a crit- 
ical phase to a columnar state with long-range order. In 
this paper we confirm that this is indeed the case. At 
finite hole density, hitherto available results arc limited 
to the form of the associated RK QDM Hamiltonia n 17 i 18 
and numerical results for small systems. 

In this work we employ analytic method s 19 ! 20 ' 21 ! 22 
combined with advanced classical Monte Carlo (MC) 
simulation s 23 ! 24 to probe the correlations in the doped 
RK wave functions, investigate the phase diagram and its 
phase transitions. The methods used here can be readily 
generalized to the case of non-bipartite lattices, for which 
a number of important results have been published.— ' 25 ! 26 
In Section mi we describe the construction of two gener- 
alized quantum dimer RK Hamiltonians that we used 
in our study. A similar but independent construction 
has been presented by Castelnovo and coworkers^ 7 - and 
by Poilblanc and coworkers 1 ^. The RK wave functions 



of these generalizations of the quantum dimer model 
have either a fixed number of holes or a variable num- 
ber of holes and a fixed hole fugacity. The ground-state 
wave functions of both models at their associated RK 
points correspond to a canonical dimer-monomer system 
in the canonical and grand-canonical ensembles respec- 
tively. Near the end of the paper, in Section IVII wc 
introduce a third Hamiltonian, with an associated RK 
wave function, to study the effects of hole interactions 
which compete with phase separation at the first-order 
transitions that we find for both models. 

In Sections IIIII and IIVI wc study the correlations and 
the phase diagram for the ground states encoded in these 
wave functions by means of an analysis of the equivalent 
classical statistical system of dimers and holes for the 
RK Hamiltonians of Section [Til hi Section IIIII we sum- 
marize the results of a mean-field theory for both non- 
interacting and interacting classical dimer models at fi- 
nite doping. The details of the mean-field theory are pre- 
sented in Appendix [Al where we compute the hole-hole 
correlation function and derive a qualitative phase dia- 
gram as a function of hole density and dimer interaction 
parameters. The main result of this simple mean-field 
theory is that the phase diagram at finite hole density 
contains two phases, a dimer-hole fluid and a columnar 
dimer solid. The columnar-liquid transition is continu- 
ous at weak coupling and turns first order at a tricriti- 
cal point. Naturally, the critical and tricritical behavior 
are not correctly described by the mean-field theory, al- 
though the general topology of the phase diagram is cor- 
rect and, remarkably, even the location of the tricritical 
point is consistent with what we find in the MC simula- 
tions of Section [Vj 

In Section ITVl we present a detailed analytic theory of 
the critical behavior of interacting classical dimers. Sec- 
tions |IV3 and |IVB] focus on the field-theoretic Coulomb- 
gas approach for this model at zero and finite doping, 
respectively. We show that, up to a critical value of a 
parameter, the undoped RK wave function describes a 
critical system with continuously varying critical expo- 
nents, with a phase transition (belonging to the 2D KT 
universality class) to a state with long-range columnar 
order. At finite hole doping we find a hole-dimer liquid 
phase (with short-range correlations) and a stable phase 
with long-range columnar order. At low hole densities 
the phase boundary is a line of fixed points with vary- 
ing exponents ending at a KT-type multicritical point 
where the transition becomes first order. We present a 
field-theoretical treatment of this tricritical point and a 
theory of the evolution of the behavior of the columnar 
and oricntational order parameters and of their suscep- 
tibilities along the phase boundary. Past the tricritical 
point the system is found to exhibit a strong tendency 
to phase separation, which we verify in our numerical 
simulations (Section |Vj. In Section [VTJ we consider the 
effects of direct hole-hole interactions near the first-order 
phase boundary, and discuss one of the many inhomogc- 
neous phases which arise in this regime instead of phase 



3 



separation. 

In Section [Vj we confirm our analytic predictions via 
extensive classical MC simulations of the generalized RK 
wave functions. For the study of the line of critical points 
at low doping we employ the canonical generalized geo- 
metric cluster algorithm (GGCA), whereas the first-order 
transition is studied via grand-canonical Metropolis-type 
simulations. The GGCA algorithm enables us to study 
relatively large systems, up to 400 x 400, for a range of 
dopings, < x < 0.06, and to investigate the finite-size 
scaling behavior. The accessible range of system sizes 
should be compared to what can be reached for full quan- 
tum models, away from the RK condition, where avail- 
able methods, such as exact diagonalization and Green 
function Monte Carlo, allow the study of only very small 
systems with few holes.— ^ 

In the undoped case we confirm the existence of 
a Kosterlitz-Thouless transition from a line of critical 
points to an ordered columnar state, as found in the work 
of Alet and coworkers . 14 ' 16 We study the scaling behavior 
of the columnar and oricntational order parameters and 
of their susceptibilities. We also use a mapping of the 
orientational order parameter of the interacting classical 
dimer model to the staggered polarization operator ob- 
tained by Baxter for the six-vertex model^i to fit our MC 
data and find an accurate estimate of the KT transition 
coupling in the undoped case. 

At low doping, we study the transition from the dimcr- 
hole fluid phase to the columnar state. We confirm that 
the scaling dimension of the columnar order parameter 
operator is equal to 1/8, as predicted by our analytic re- 
sults of Section llVBI We also present a typical set of data 
that demonstrates how the scaling dimension of the orien- 
tational order parameter varies, again in agreement with 
the analytic results of Section llVBl and use these results 
to locate numerically the phase boundary. We then turn 
to the behavior at larger doping and stronger couplings 
where the transition becomes first order. We study this 
regime using grand-canonical Metropolis Monte Carlo 
simulations. We confirm the first-order nature of the 
phase transition by means of a careful analysis of the 
finite-size scaling behavior of the order parameters across 
the phase boundary and of their susceptibilities. We use 
these results to locate the phase boundary in the first- 
order regime as well. In Section IVI[ we use MC simula- 
tions to study the effects of a direct hole-hole repulsion 
which suppresses the effects of phase separation, leading 
instead to a complex phase diagram of inhomogeneous 
phases, of which we only study its most commensurate 
case. 

In Section fVIIl wc study the elementary quantum ex- 
citations of the doped QDMs satisfying the RK condition 
using the single-mode approximation. We only present 
the main results and have relegated the details to Ap- 
pendix [Bj We find that in the dimer-hole liquid phase, 
hole and dimer density fluctuations have quadratic dis- 
persions E(k) ~ k 2 . Thus this phase should be charac- 
terized as a Bose- Einstein condensate of bosonic charged 



particles (holes), but not really a superfluid, for reasons 
similar to those of Rokhsar and Kivelsonpi We summarize 
our overall conclusions in Section [Villi 

While this manuscript was being completed (and ref- 
ereed) a number of independent studies of aspects of this 
problem have been published ] 16 ' 17 ' 18 Our results agree 
with those in these references wherever they overlap, as 
noted throughout this paper. 

II. QUANTUM HAMILTONIANS FOR 
INTERACTING DIMERS AT FINITE HOLE 
DENSITY 

The Hamiltonian of the quantum dimer model (QDM) 
can be written in the Rokhsar-Kivelson (RK) formi as 
the sum of a set of mutually non-commuting projection 
operators Q p , 

H = Y^Qp, (2-1) 

M 

where {p} denotes the set of all plaquettes of the square 
lattice. Each projection operator Q p acts on the states of 
the dimers and holes of a plaquette p (or set of plaquettes 
surrounding p) . In the simplest casei each Q p acts only 
on the states labeled by the dimer occupation numbers 
of the links of the plaquette p. In this case, the ground 
state is described by the short-range RVB wave functions*. 

|RVB)=]T|C), (2.2) 

{C} 

where {C} is the set of (fully packed) dimer coverings 
of the 2D square lattice, and {\C)} is a complete set 
of orthonormal states. If one regards the dimers as 
spin singlet states (with the spins residing on the lat- 
tice sites) each configuration represents a set of spin 
singlets or valence-bonds The dimer representation 
ignores the over-completeness of the valence bond sin- 
glet statesi 1 - This problem can be made parametrically 
small using a number of schemes, including large N 
approximations^ and decorated generalizations of the 
spin-1/2 Hamiltoniansi^ 

It is possible and straightforward to generalize the 
QDM construction so as to include other types of in- 
teractions and coverings. In Ref. [ll] it was shown how 
to extend this structure to smoothly interpolate between 
the square and the triangular lattices. It was also shown 
there that the same ideas can be used to construct a 
quantum generalization of the two-dimensional classical 
Baxter (or eight-vertex) model. In all of these cases, 
the RK form of this generalized quantum dimer model 
has an exact ground-state wave function whose ampli- 
tudes are equal to the statistical (Gibbs) weights of an 
associated two-dimensional classical statistical mechan- 
ical system on the same lattice. Thus, if the classical 
problem happens to be a classical critical system, the 
associated wave function now describes a 2D problem 
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FIG. 1: (color online) Illustration to clarify the construction 
of the Hamiltonian (|2.3p . The central dimer pair is present 
in all configurations d and C'i, whereas the four surrounding 
dimers may or may not be present. The index i enumerates 
the 2 4 possible arrangements of the surrounding dimers. The 
configuration C' differs from the corresponding C by a single 
flip resonance of the central dimer pair. 



at a quantum critical point. In Ref. [XT| such quantum 
critical points were dubbed "conformal quantum critical 
points" since the long-distance structure of their ground- 
state wave junctions exhibits 2D conformal invariance. 
Here we are interested in a different generalization of the 
QDM in which wc consider dimer coverings (although 
not necessarily fully packed) of the square lattice. We 
will also consider 2D Hamiltonians whose wave functions 
correspond to classical interacting 2D dimer problems 
with local weights. Similar but independent construc- 
tions have also been proposedi 17 ' 18 i 30 

Trying to be as physical and local as possible, we keep 
the quantum-resonance terms as simple as before (single 
plaquettc moves), but the potential terms (which again 
have a central plaquette^) now have fine-tuned couplings 
that depend on the nearby plaquettes. Explicitly, we 
have (cf. Fig.rjJ: 



c,)(c ? '|-|c:)(a 



w 



-Rc, 



(2.3) 



where Rc i and Rc. denote the number of pairs of present 
dimers formed in configurations d and C[ respectively 
The Hamiltonian of Eq. (|2.3|) is designed in such 
a way that it annihilates any superposition of dimer- 
configuration states which have amplitudes that are of 
the form w Np , where w is the parameter appearing in 
the Hamiltonian and N p is the number of pairs of neigh- 
boring dimers in the configuration*^ In this sense, the 
Hamiltonian is a sum of projection operators and conse- 
quently there is a unique ground state for each topologi- 
cal sector which must be composed of the superposition 
of these specially weighted configurations. The ground- 



state wave function, \G), the state annihilated by all the 
projection operators, for this system is: 



|G) = 



1 



where Z(w 2 ), the normalization of this state, 
Z(u?) = J2™ 2Np[C] 

{C} 



(2.4) 



(2.5) 



has the form the partition function of classically interact- 
ing dimers with a coupling u = — 21nu> between parallel 
neighboring dimers. In the following, wc will assume an 
attractive coupling, u < or w > 1. The case u > was 
studied for the fully-packed case in Ref. QjJ 

There are two different ways in which we can add dop- 
ing to our system, while still being able to determine the 
ground state. If we add the following fine-tuned hole- 
related terms to the initial Hamiltonian 



H. 



hole 

canonical 



"£hole ^2 \ 
<ijk> I 



. !)<_• 



h.c. 



. IX. I 



(2.6) 



then the resulting ground state becomes 



\G 



1 



N h 



J2 w N P^\\C Nh ), (2.7) 



where the number of holes Nh is now fixed at a specified 
value. The norm of this wave- function, Z(w 2 ), is the 
canonical partition function for the set of dimer coverings 
with a fixed number of holes. 

On the other hand, if we add the following terms, which 
do not conserve the number of holes in the system, to the 
Hamiltonian (|2.3p 



rrholc 

grand— canonical 



"iholc ^ < 
links y 



(2.8) 



then the ground state of the system becomes 



\G)= . * 2 T W C 'lO. (2-9) 
y/Z(w 2 ,z 2 ) £rf 



with 



2N p [C]2N h [C] 



z(w 2 ,z 2 ) = u>- 

{C} 



(2.10) 



Equation (|2.8|) includes a natural, short-range repul- 
sion between neighboring holes and an off-diagonal term 
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which represents creation-annihilation of dimers. Fur- 
thermore, there is a dimer fugacity term, as in every per- 
turbativc derivation of a quantum dimer model.— We note 
that Eq. (|2.10[) has the same form as a grand partition 
function for dimer coverings of the square lattice. This 
partition function now depends not only on the inter- 
action u defined below Eq. (|2.5p . but also on the hole 
chemical potential (i/\u\ = 21nz. 

Since the canonical and grand-canonical ensembles be- 
come equivalent in the thermodynamic limit, the two 
ground-state wave functions (|2.7|) and (|2.9p must corre- 
spond to the same ground-state physics. Furthermore, it 
is clear that for w = 1 the models are located at the usual 
RK point of the quantum dimer model on the square lat- 
tice. For a system with periodic boundary conditions, 
each configuration C contains only an even number of 
holes, with half of the holes on either sublatticc. 

We remark that the fact that the exact ground-state 
wave function is a sum (as opposed to a product) of the 
ground states of sectors labeled by the number of holes 
on the lattice, is due to the resonance term that we have 
used to represent the motion of holes. In particular, we 
have assumed that a dimer can break into two holes which 
themselves repel each other. In the limit of very strong 
hole-hole repulsions, in strong-coupling perturbation the- 
ory, it is straightforward to recover a fixed-hole density 
sector with a single-hole resonance move involving three 
sites in any direction; the coupling strength in Eq. (|2.6[) 
then becomes ihoic ~ z 4 ihoic, and thus, in the limit z — > 0, 
it reduces to an effective hopping amplitude for the holes. 



III. MEAN-FIELD RESULTS 

To examine the physics described by the ground-state 
wave functions obtained in the previous section we start 
with a discussion of a mean-field theory of the phase 
diagram. We use the standard approach of regard- 
ing the probability densities, obtained by squaring the 
wave function, as the Gibbs weights of a classical two- 
dimensional system and focus on the interacting dimer 
model on the square lattice at finite hole density. Al- 
though mean-field theory is insufficient to describe two- 
dimensional critical systems, it is a useful tool to obtain 
qualitative features of the phase diagram as well as the 
behavior deep in the phases, away from criticality. 

The details of this theory are presented in Appendix lAl 
We begin by constructing a Grassmann representation of 
the partition function for an interacting dimer model us- 
ing the standard methods introduced by SamueL 31 i 32 The 
resulting theory involves Grassmann (anti-commuting) 
variables residing on the sites of the square lattice. The 
action of the Grassmann integral is non-trivial and is 
parametrized by the hole fugacity z and a coupling be- 
tween dimers V = z 2 (e u — 1), where u = — 21nu>. 

Since the action of the resulting Grassmann path in- 
tegral is not quadratic in the Grassmann variables, it 
cannot be reduced to the computation of a determinant. 



Thus, we use a standard mean-field approach which, 
in this case, involves the introduction of two Hubbard- 
Stratonovich (bosonic) fields fa and Xij> defined on the 
sites and links of the square lattice respectively. Upon 
integrating out the Grassmann variables one obtains an 
effective theory for the fields fa and Xij which, as usual, 
is solved within a saddle-point expansion. The dimer mo 
and hole n densities, as well as the columnar order pa- 
rameter m, can be expressed straightforwardly in terms 
of the fields fa and Xij ■ From this effective theory one can 
compute an effective potential T and the configurations 
of the observables of interest, n, mo and m, as functions 
of z and V, and determine the phase diagram. 

As a function of hole density (or hole fugacity) and u, 
we find that the phase diagram has two phases (shown 
qualitatively in Fig. [21) , namely a dimer- hole liquid phase 
and a hole-poor phase with long-range columnar order. 
The nature of the transitions between these phases is in- 
correctly described by the mean-field theory, particularly 
at zero doping and near the tricritical point. The cor- 
rect behavior is the subject of a detailed analysis in the 
subsequent sections. Nevertheless, the mean-field phase 
diagram correctly predicts that at low hole densities and 
moderate values of u the transition between the dimer- 
holc liquid and the columnar solid phase is continuous, 
that for large u the transition is first order, and that there 
is a tricritical point at 

utr ^ 2.733, z tr ~ 0.075 . (3.1) 

Remarkably, the Monte Carlo simulations presented in 
Section [V] yield a tricritical point at a location quite con- 
sistent with these values. 

Another correct prediction of the mean-field theory 
is the behavior of the connected hole density correla- 
tion function deep in the dimer-holc liquid phase. This 
prediction, also discussed in detail in Appendix [A] fits 
surprisingly well the Monte Carlo simulations of Krauth 
and Moessner^ 3 - performed for a system of classical non- 
interacting (u = 0) dimers at finite hole density. The 
mean-field result is consistent with the simulations for a 
quite broad region of densities even quite close to the fully 
packed limit z — > where, naturally, a wrong correlation 
length exponent is predicted. 

IV. PHASE DIAGRAM AND CORRELATIONS 
FOR INTERACTING DIMERS AND HOLES 

We now turn to a more precise analysis of the phase 
transitions of the interacting classical dimer models as 
a function of hole density on interaction parameter u. 
Here we take advantage of a wealth of information and 
methods from two-dimensional systems, exact solutions 
and conformal field theory, to analyse the behavior in 
detail and extract conclusions that will be quite useful 
for the analysis of the wave functions. We begin with 
a discussion of the undoped case, and then discuss the 
physics at finite hole density. 
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FIG. 2: Qualitative phase diagrams for the interacting doped 
dimer model, a) Phase diagram as a function of the fugacity 
z for the holes and u = — In w, the dimer interaction coupling 
constant. The line of continuously varying exponents for low 
fugacity evolves to a line of first-order transitions, b) Phase 
diagram in terms of the hole density ph and the coupling 
constant u. Here the line of first-order transitions opens up 
into a two-phase coexistence region. 



A. Interacting Dimers at Zero Hole Density 

It is a well known fact that both classical and quan- 
tum two-dimensional dimer models can be represented in 
terms of height models. For the classical case, this map- 
ping is well know n n i 33 ' 34 i 35 . The mapping for the quan- 
tum case has also been discussed extensively^ ' 11 ' 12 ' 36 ' 37 . 
In both cases the mapping relates each dimer configura- 
tion to a configuration of a set of integer-valued (height) 
variables, h(r), which reside on the sites of the dual lat- 
tice, a square lattice in the case of interest here. Thus, 
this mapping amounts to a duality transformation. 

An alternative picture follows from realizing that 
dimer configurations can be mapped onto the degener- 
ate ground-state configurations of fully frustrated Ising 
modcl o 33 ' 38 . In our case, the corresponding spin model 
is the fully frustrated Ising model on the square lattice 
(FFSI) at zero temperature. In the Ising model picture, 



each dimer is dual to an unsatisfied bond of the fully 
frustrated Ising model and classical dimer interactions 
correspond to second neighbor interactions in the square 
lattice FFSI 3 ^. It is easy to see that holes correspond to 
unfrustratcd plaquettes in the FFSI. In this work we use 
primarily the language of the height representation. 

To map the square-lattice classical dimer model onto a 
height model, we follow the prescription given in Ref. [Til 
One first assigns a height variable to each plaquctte. In 
going around a vertex on the even sub-lattice clockwise, 
the height changes by +3 if a dimer is present on the 
link between the plaquettes, and by —1 if no dimer is 
present on that link. On the odd sub-lattice, the heights 
change by —3 and +1 respectively. The dimer constraint, 
that one lattice site belongs to one and only one dimer, 
implies (for the square lattice) that, for the fully packed 
case, there are four possible configurations of dimers for 
each lattice site. In the dual height model this is repro- 
duced by the period four property h = h + 4 of the al- 
lowed height configurations. It is easy to see that for the 
allowed configuration the average values that the height 
variables can take at a given site of the direct lattice (a 
vertex) are ±3/2, ±1/2. On the other hand, a uniform 
shift of all the heights by one unit, h — > h + 1, leads to 
an equivalent state. This mapping works strictly speak- 
ing only for the fully packed case. Holes are sites that 
don't belong to any dimer and thus represent violations 
of the full packing rule. They play the role of topological 
defects ("vortices") in the (dual) height representation. 

The exact solution of the non-interacting fully packed 
dimer model on the square (and other) lattice has 
been known for a long time^. In particular the long- 
distance behavior of the dimer density correlation func- 
tions and the hole density correlation functions arc 
known explicitl y 40 ' 41 . These correlation functions obey 
power law behaviors and show that this is a critical sys- 
tem. Here we will use the standard approach to map the 
exact long distance behavior of two-dimensional critical 
systems to the behavior of the simplest critical system, 
the Gaussian or free boson models. This approach is 
consistent for the free dimer model on an even lattice with 
periodic boundary conditions since its central charge (or 
conformal anomaly) is also c = 1. 

We will use as reference states ("ideal states" in the 
terminology of Kondev and Henley^ 2 -) the four colum- 
nar states of the dimer coverings, which have the largest 
number of flippable plaquettes, and use them to define 
an effective field theory for this problem^. We will as- 
sign a uniform value to the coarse-grained height field 
h = 0,1,2,3 to each reference (columnar) state. Let 
n x {T) and n y (y) represent the coarse grained dimer den- 
sities of the horizontal link with endpoints at the pairs of 
lattice sites r and f+e x , and vertical links with endpoints 
r and r + e y respectively. Here e x and e y are two lattice 
unit vectors along the x and y directions respectively, 
with a lattice spacing of 1. 

We can now define the columnar local order parameter 
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as the two component vector 

O x {x,y) = n x (x,y) -n x (x + l,y) 

O y (x,y) = n y (x,y)-n y (x,y + l) (4.1) 

which clearly correspond to the (7r, 0) and (0, 7r) Fourier 
components of the dimer densities. This two-component 
order parameter takes four distinct values for each one 
of the columnar states, and changes sign under shifts by 
one lattice spacing in either direction. It also transforms 
as a vector under 90° rotations. Thus it is the order 
parameter for columnar order. 

1. Effective Field Theory: the non-interacting case 

The fluctuations of the free field h(r) are described by 
a continuum Gaussian (free boson) model. We will find 
it simpler to work with the rescaled height field <fi = ^h. 
For this field the periodicity condition h — > h + 4 be- 
comes <j> — > <fi + 2ir. (For the rescaled field the ideal states 
are <f> = 0,tt/2,tt,3tt/2.) Thus, the allowed operators 
are 2tt periodic functions of <p, and are either derivatives 
of 0, or the exponential (or charge) operators exp(±i(/>), 
exp(±2i(/>), exp(3«0) and exp(±4i0), which are 2tt peri- 
odic functions of <j>. 

The action S for the rescaled field is 

S = J d 2 *| (V0) 2 (4.2) 

For the free dimer model the stiffness is K = -r- (see 
below) . 

By matching the exact correlation functions of the 
free dimer model on the square lattice one readily finds 
the following operator identification of the coarse-grained 
dimer densities in terms of free field operators^ 

n x -\ = J-(-l)-+^+i[(-ire^ + c.c] (4.3) 

4 Z7T Z 

n y -\ = J-(-l)^+^+i[(-irie^ + c.c] 

4 Z7T Z 

(4.4) 

In Ref. d it was shown that this is an operator iden- 
tity for the free dimer model on the square lattice in the 
sense the the asymptotic long-distance behavior of the 
dimer density correlation functions computed with this 
Gaussian model are the same as the exact long distance 
correlation functions for the free dimer problem on the 
square lattic e 40 ' 41 provided the stiffness K = j-. Notice 
that, with this identification, when the field <fi takes each 
of the values 0, tt/2, it, Ztt/2 (the "ideal states") the den- 
sity operators take four distinct values which reflect the 
broken symmetries of the four columnar states. 

From the operator identification of Eq. (|4.4|) the colum- 
nar order parameter is, up to a normalization constant, 

02, = cos 0, y = sin0 (4.5) 



Due to the effects of dimer-dimer interactions the form 
of this effective action is 

S = J tfx y (V0) 2 + perturbations (4.6) 

The effect of the interactions is a finite rcnormalization of 
the stiffness K away from its free dimer value, K = 

We saw above that, due to the dimer constraints, the 
allowed charge operators arc 0„(f) = cxp(m0(r)). We 
also saw that the columnar order parameter is propor- 
tional to the operator ^(Oi(f) +0_i)(r) = cos </>(?), and 
carries the unit of charge n = 1. One can also define 
vortex or magnetic operators^, and example of which is 
the hole. A vortex operator causes the field </> to wind 
by 27rm, where m is the vorticity (or magnetic charge). 
One can similarly define a general composite operator 
O n .m{r) with n units of (electric) charge and m units of 
vorticity (or magnetic charge). Its scaling dimensions, 
A(n, m) arei^ 

A n , m (K) = -P--+TrKm 2 (4.7) 

We can now use these results to identify a few opera- 
tors of interest and give their scaling dimensions. These 
results are summarized in Table U 

1. The columnar order parameter is the elementary 
charge operator O±i,o an d has no vorticity. On the 
(columnar) ideal states 0, ir/2, ir, and 37r/2 this 
operator takes the values 1 ,i ,— 1, and —i respec- 
tively. Its scaling dimension is Ai^-ftT) = At 
the free dimer point, K = j- , its scaling dimen- 
sion is Ai i o(j-) = 1. This is consistent with the 
exact result o 40 ' 41 that the density correlation func- 
tion falls off as 1/r 2 . The operator identification of 
Eq. (|4.4[) is based on these facts^. 

2. The operator O±2,o = exp(±2i0) takes the values 
1, — 1, 1 and —1 on each of the ideal columnar 
states. It is clearly the order parameter for sym- 
metry breaking by 90° rotations: it is the order 
parameter for orientational symmetry. 

3. The operator with the highest allowed electric 
charge is O±4,o = exp(±4£<^>). Its scaling dimension 
is A^ t o(K) = At the free dimer point it has 
dimension A4 i o(t— ) = 16, and it is a strongly irrele- 
vant operator. This operator arises naturally due to 
the fact that the microscopic heights h take integer 
values, and hence height configurations which differ 
by an uniform integer shift are physically equiva- 
lent. This operator does not break any physical 
symmetry of the dimer model. 

4. The hole operator is represented by the fundamen- 
tal vortex operator Oo,±i- A vortex with unit pos- 
itive magnetic charge corresponds to a hole on one 
sublattice, and a vortex with unit negative mag- 
netic charge to a hole on the other sublattice. The 
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scaling dimension of the vortex (hole) operator is 
A 0j i(A") = ttK. At the free dimer value, the scal- 
ing dimension of the hole operator is A^id^r) = \, 
which is consistent with the exact result that the 
hole-hole correlation function decays at large 

distances^. 

5. The operator which describes a pair of holes on 
nearby sites of the same sublattice is represented 
by the operators Oo,±2 which carry two units of 
magnetic charge (vorticity). This operator cre- 
ates (or destroys) a diagonal dimer connecting 
nearby points on the same sublattice^ In the 2 + 1- 
dimensional quantum dimer model, this operator is 
useful to describe the possible pairing of holes. This 
operator has dimension Ao,2(A) = AttK. Its scal- 
ing dimension at the free dimer point is Ao.2(^r) = 

1, and is relevant for A < 77-. As noted in Refs. |251 
I — I — 
and[ll|, this operator maps the square lattice into a 

deformed triangular lattice. The irrelevancy of this 

operator for A > ^ implies that this line of fixed 

points also exists for a deformed triangular lattice, 

as discussed recently in Ref. [26l 

6. The free dimer problem is a free fermion system, 
which can be solved by Pfaffian methods i 25 ' 31 i 40 
This is actually a theory with two free real (Ma- 
jorana) fcrmions or, equivalcntly, one free com- 
plex (Dirac) fermion, at its (massless) fixed point, 
whose central charge is also c = 1. The appro- 
priate fermion operator is a composite operator of 
the order and disorder operator o 19 ' 43 which, in this 
case, is Oi v At the free dimer point the fermion 

operator has scaling dimension Ai^(^) = \ and 

(conformal) spin nm = i (as it should for a free 
fermion). At particular values of K it is also possi- 
ble to define parafcrmion operator o 44 ' 45 operators 
which obey fractional statistics. For instance, at 
the KT point, K = -| (see below), the opera- 
tor 1 1 has dimension j and spin j, and it is a 
semion. In fact, and not surprisingly, a fermionic 
approach^ 5 - can be used to map this critical line 
onto an Euclidean version of the Luttinger-Thirring 
model. The coarse-grained height model descrip- 
tion we sue here corresponds to the bosonization 
approach of the fermionic version of this problem. 



2. Effective Field Theory: interactions and phase transition 

We now turn to the effects of dimer-dimer interactions. 
We recall that we are considering only interactions of a 
pair of dimers in a plaquette. The interaction energy is 



^ \n x (f^n x (r + e y )+n y (r)n x (f+e x ) (4.8) 



We wish to find the corresponding operator in terms of 
the coarse grained (rescaled) height field <fi. This can 





columnar 
Oi.o 


rotationa 
C»2,o 


C»4,o 


hole 

Oo,i 


hole pair 

O ,2 




1 


4 


16 


1 

4 


1 


TT 


1 

8 


1 
2 


2 


2 


8 



TABLE I: Scaling dimensions of the order parameters 
(charge) and hole (vortex) operators at the free dimer point, 
K = j-, and at the KT point, K = -. 



be done by using the operator product expansion (OPE) 
of the coarse-grained form of the dimer density opera- 
tors, Eq. (|4.4p . in terms of the field cf>. We will also 
need the standard OPE of the (normal ordered) charge 
operator o 21 ' 46 ! 47 



: cos(n<^>(a;)) : : cos(n(f>(y)) 
1 



cos(2n</>(x)) 



(^|s-y| 2 )» a 

sin(n0(x)) 
1 

2\n 2 



l-n 2 \x-y\ 2 : (V<£) : +• 



(4.9) 



sin(n</>(y)) := — — : cos(2n</>(:r)) 



(jj?\x-y\ 2 ) 



l~n 2 \x-y\ 2 : (V0) : +. 



(4.10) 



where fi is a short-distance cutoff and the ellipsis rep- 
resents the contributions of irrelevant operators. Using 
these results we find that the net effect of the interactions 
is to rcnormalizc the stiffness A upwards 



K 



i_ + i( 1 + _L]„ + cV) 



(4.11) 



where we have denoted u > for attractive interactions. 
This expression for the renormalization of A is only ac- 
curate to linear order in the dimer-dimer interaction. 
Higher order renormalizations (in u) would result if the 
effects of irrelevant operators are also taken into account. 
The relation between K and the microscopic model is 
non-universal and can only be determined either from 
an exact solution or from a numerical simulation. One 
can determine the function A (it) from the Monte Carlo 
simulations we present elsewhere in this paper. What is 
important is that these non-universal effects affect only 
the relation between the coefficients of the effective the- 
ory and not the form of the effective theory itself. Thus, 
the effective action of the field theory for the interacting 
classical dimer model at zero hole density has the form 



S 



d 2 x 



y (V0) 2 + ff COs(40) 



(4.12) 



where we have included the effects of the charge 4 per- 
turbation, cos(4</>) = cos(27iA), which biases the coarse- 
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grained height field to take integer values. 

For an anisotropic dimcr-dimcr interaction, which 
arises form a term which weights differently the inter- 
actions between parallel horizontal dimers from those 
of parallel vertical dimers, we would have also found a 
cos(2<^>) operator in addition to an anisotropy for the stiff- 
ness. Thus, an anisotropy in the dimer-dimer interaction 
is a relevant perturbation which couples to the orienta- 
tional order parameter cos(2<^). 

In Table |T] we see that as the attractive interactions 
grow there will be a critical value of the interaction u at 
which the stiffness K(u c ) = ^. At this point, the cos(40) 
operator has scaling dimension A 2i o = 2, where it be- 
comes marginal. For u > u c (K > K(u c )) this operator 
becomes relevant. We also see that at K{u c ) = — the 
columnar order parameter, cos (f>, has scaling dimension 
1/8 and it is the most relevant operator in this problem. 
Thus this is a phase transition from a critical phase, for 
K < — , without long-range order but with power law 
correlations, to a phase with long-range columnar order, 
for K > — , in which the columnar order parameter has 
a non-vanishing expectation value. This is a standard 
Kosterlitz-Thouless (KT) transitioni&i^ which is nat- 
urally described by the sine-Gordon field theor y 51 ' 2 ' 53 ' 
whose (Euclidean) action is given in Eq. (|4.12p . The 
only difference between this problem and the standard 
KT transition, e.g. the classical 2D XY model, is that 
the phase with a finite correlation length is ordered: it is 
a columnar state with a four-fold degenerate non-uniform 
state, whereas the finite-correlation length phase of the 
XY model (and of its dual surface roughening model) has 
a non-degenerate translation invariant state. In spite of 
these global differences, this phase transition is in the 
KT universality class. Thus the well known behavior of 
the correlation functions at the KT transition apply to 
this case as wel l 14 ' 16 ' 17 ' 18 . In Section [V] we verify this be- 
havior by a detailed Monte Carlo study of the columnar 
and oricntational order parameters and of their associ- 
ated susceptibilities for the interacting dimer model. 



magnetic charges in two dimension o 55 ' 56 



Z(w, z) = ex P 

{n(r),m(R)} 



N 2 
AttK 



n{f) In \f — r '| n(r') 



-nK m (R) hi\R-R'\ m(R') 



x exp 



Y^lnw n(f) 2 + ^lnz m(R) 2 



:J2N n(r) m{R) 9(r - R) 



r,R 



(4.13) 



where prime denotes that the sum is restricted to neutral 
configurations with vanishing total charge and vanishing 
total vorticity, i.e. J2r n (^) = Z^R m (R) = 0- Here 
0(r — R) is the angle between a vortex at R, as seen from 
a charge at r, measured with respect to the (arbitrary) 
x axis. It is the Cauchy-Ricmann dual of the logarithm, 



G(f-r') 

e(f-f') 



— In \r — f" 

. i'y 



tan 



.r 



-V 2 G(r-r") = 2tt 5 2 {r-r') 



(4.14) 

(4.15) 

(4.16) 
(4.17) 



As usual, the logarithmic interaction is regularized so 
that it vanishes for f ' = f 1 ' . The short distance behavior 
of the interactions is absorbed in the fugacities z and w. 

For the case we are discussing here, the GCG of inter- 
est has N — 4, and Eq. (|4.13|) is just the the Coulomb-gas 
form of the partition function of the Z4 model, and for 
the related 2D Ashkin- Teller model. This is a well un- 
derstood syste m 21 i 22 ' 50 i 52 and in our case, it corresponds 
to the grand-partition function for the doped interacting 
dimer model on the square lattice at low hole densities. 

In the limit of low fugacities, z <C 1 and w <C 1, the 
GCG is equivalent to a (generalized) sine-Gordon field 
theory in two-dimensional Euclidean space-time, whose 
(Euclidean) Lagrangian is given b y 21 i 52 



B. Interacting Dimers at Finite Hole Density 



We now consider the dimer model at finite hole den- 
sity, away from the full packing condition. The classical 
partition function for this problem is given in Eq. (12.101) , 
where the weights (fugacities) z and w are related to 
the coupling constant u and the chemical potential fi 
as described earlier. Recall that the 2D classical par- 
tition function Z(w, z) is the norm of the ground-state 
wave function \G), of Eq. (|2.4p . of the 2D doped quan- 
tum dimer model. In terms of a sum over configurations 
of electric charges n and magnetic charges (vortices) to, 
the partition function Z(w , z) is equivalent to that of a 
generalized (neutral) Coulomb gas (GCG) of electric and 



C E = - {d^f 



Zz / iv \ Zw / r— 

— cos — = 4>\ — cos I 27TV K 

a 2 \y/K J a 2 \ 

(4.18) 

Here </> is the dual field, defined by 



e^dvef) = id^ (4.19) 

In Eq. (|4. 18|) we can see by inspection that the oper- 
ators cos(^= 4>) and cos{2-k^/K 4>) can be identified 
respectively with the operators ^(On.q + O-n.o) an d 
5(00,1 + Oo, -1) discussed above. 

We will use this effective field theory to study the tran- 
sition between the liquid and the ordered phases of the 
interacting dimer model. At z = w = this is the KT 



10 



transition discussed above. For general N, both oper- 
ators have the scaling dimension if -Jjjl— = 2-k\J~K^ i.e. 

the theory is self-dual, which happens for K = For 
N = 4, K = — , both operators have scaling dimension 
2 and both are marginal. This is the only case we will 
discuss here. (A detailed discussion of the more general 
case of N > 4 was given by Lecheminant et al*£) 
For N = 4 the Euclidean Lagrangian becomes^ 2 - 



i(t + x) and z — r — ix = i(t — x), in imaginary and real 
time respectively. 

It is easy to se o 20 ' 22 ^ 57 that the following dimension 1 
chiral operators 



J't 



T~ 



'2tt 



d z (j)L 
dz4>R 



J 



J 



± _ 1 . e ±iV8ii(j>L 



R 



2ir 



(4.26) 



1 2 



2z 



■ cos 



'87T (j) ) cos I 



2w 
~ 2~ 



a 



(V8tt (j>) 
(4.20) 

It turns out that this a problem which can be solved 
exactl y 20 i 22 The most direct way of doing this is to per- 
form an analytic continuation from 2D Euclidean space 
to 1 + 1-dimensional Minkowski space-time, i.e. to think 
of this problem as a 1 + 1-dimensional quantum field the- 
ory The Hamiltonian density of the equivalent 1 + 1- 
dimcnsional field theory is 



with Jf; R = R ±iJ^ R , are the generators of an su(2)j 
Kac-Moody algebra given by the OP E 22 ' 57 



r L (z)j b L {w) 



Sab 



€abc 



8it 2 (z-w) 2 8tt 2 (z-u;) 

&ab ■ ^-abc 

8ir 2 (z-w) 2 8ir 2 (z-w) 



Jr.(u>) 



(4.27) 



where a,b, c = x,y,z and e a b c is the Levi-Civita tensor. 
We also note for future use the following operator iden- 



— ^ cos (V8ir <p*j — 



2w 



-2 
+2 



(z + w) 

2 C0 

or 

(z — w) . 







tifications (still for K =2/%) 




cos ^ 




e i\f&K(j) 


= o 4 . 

= O ,i 


cos(-> 








sin(\ 


^l) (4.21) 




= T2,±i 



(4.28) 



where we have used the fact that II, the canonical mo- 
mentum conjugate to the field <p is simply related to the 
dual field 6, 



The Hamiltonian of Eq. (|4.21[) can be written in terms 
of these generators and has the Sugawara form 5 - 7 -: 

2tt 



n = ~~FT~ = d t<t> = ~&x 

0(f) 



(4.22) 



and obey equal-time canonical commutation relations 

[4>{x),n{y)]=i8{x-y) (4.23) 



In Eq. (|4.2ip we used the decomposition of the field (f> 
and the dual field <j> into right and left moving fields (f>n 
and 4> L , 



4>l = -OH 

whose propagators are 



= <PL ~ <PR 

^R=\{<t>-4>) (4.24) 



{4>l{z)4il{w)) = — — hx(z - w) 

47T 

(Mz)Mw)) = ~ - w) (4.25) 

47T 

where we have used the complex coordinates (not to be 
confused with the coupling constants!) z = t + ix = 



T~t = — \^Jl ■ Jl + Jr. ■ Jr. 

-8tt 2 (z + w)J x L J x R - 8tt 2 (z - w)J v L J v R (4.29) 

Thus, one recovers the old resul t 20 ' 51 that at the KT 
transition the system has an effective (dynamical) SU(2) 
symmetry. 

Given the existence of an su(2)\ symmetry one expects 
to find, in addition to the "spin one" (vector) representa- 
tion (the su(2) currents), operators associated with the 
spin 1/2 representation of sm(2)i. In the Wess-Zumino- 
Witten (WZW) version of this theor y 57 ' 58 , there is an 
operator with this property, the field g(z,z) of the WZW 
model. This is a 2 x 2 matrix-valued field with scaling di- 
mension 1/2. In our theory we also have an operator with 
scaling dimension 1/2, the operator 02,0 ~ exp(i>/2n(j)) 
which we saw was related to the order parameter for bro- 
ken rotational symmetry; see Table [I] Thus, the opera- 
tors of the spin 1/2 (spinor)representation of su(2)i are 
identified with the operators 



2tt0 



2tt0 



2tt0 



2tt0 



0-2,0 O .-i 
O 2 ,o O i 



(4.30) 
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Following the approach of Refs. I20II22I we will perform 
a global SU(2) rotation by tt/2 about the y axis which 
maps Jf >fl — * J LtR , J LtR 



after which the Hamiltonian becomes 



-Jl R and J\ R 



J L,R> 



H 2n 

T 



-8tt 2 (z + 10) J[ J Z R - 8tt 2 (z - w)J v L J v R (4.31) 

Once again, one can introduce a new bose field, which 
we will call <£>, and its dual <&, to use the representation 
of the su(2)\ current algebra of the form of Eq. (|4.26p . 
Using that z $l — ~id x $L and d s $ R — id x <& R , we can 
rewrite the Hamiltonian as 

H = Ho + Hpert 

Ho = \{d x $? + \{dj>f - 4ir(z + w)d x $ L d x $ R 



H 



pert 



W ' sm(V8n$ L ) sin(\/8~$ii) (4.32) 



Thus, along the phase boundary line z = w the term 
7i per t is absent and we see that the effective Hamiltonian 
Ho involves only marginal operators . 20 ' 22 ! 50 ' 52 

Similarly, the operators in the spin-1/2 representation 
transform as an su{2) spinor under a 7r/2 rotation about 
the y axis in su(2)\, leading to the identifications 





* O 1 = 

u, 2 


e iV2~7r$ 




O 2 ,o - 


* -°o-i = 

u, 2 


_ — i\/27T$ 




O _i - 

u, 2 


* -O 2 , = - 


e i\/27r$ 




O i - 

u, 2 


+ -0-2,0 = 


_ e -iV2^$ 


(4.33) 



where the operators on the right hand side of Eq. (|4.33[) 
are vertex operators written in terms of the new bosons 
$ and <E>. Notice that this is a duality transformation. 

However, there is an operator in this theory, namely 
O±i,0j which does not have su(2)i quantum numbers. 
As a result it will have a quite different behavior along 
the phase boundary. In contrast, all the other opera- 
tors in this theor y ca rry su(2)i quantum numbers. (As 
remarked in Ref. [2CJ this is not truly an sit(2)i theory, 
although it contains it, but rather connected with a Z2 
orbifold as well.) It is straightforward to check that the 
operators O±\ t o do not have an OPE with the opera- 
tors which do transform under su(2)i (or, rather, that 
the OPE involve only irrelevant operators.) Since the 
marginal operator which deforms the su(2)\ theory along 
the phase boundary does transform under su(2)i, the op- 
erator 0±\fl will not mix (in the sense of its OPE) with 
the marginal operator either. We will see that this im- 
plies that the dimension of this operator remains equal to 
1 /4 along the entire line, a result that is also well known 
(see for instance Ref. [53-) In contrast, the operators of 
the spin-1/2 representation, Eq. I|4.30[) . do transform un- 
der su(2)i, a fact which is generated by their OPE's with 



the su(2)i generators^ and we will now see that their 
scaling dimensions do change along the phase boundary 
line. In Section[V]we present evidence from Monte Carlo 
simulations in support of both statements. 

We can now use this approach to solve this problem ex- 
actly along the phase boundary. Formally the Hamilto- 
nian 7io of Eq. (|4.32j) . is equivalent to a spinless Luttingcr 
model with attractive backscattering interactions^. As 
in the case of the Luttingcr model, the problem is solved 
by means of a Bogoliubov transformation of the right and 
left moving bosons. This procedure breaks the su{2)\ 
symmetry explicitly. We introduce a new bose field x 
and its dual field The left and right moving compo- 
nents of these fields, \l and \Ri are linearly related to 
the left and right moving fields <f>£ and Q R by 



XL 
XR 



1 



* i+ 2 



1/1 



v/kJ $r 

1 



where k is given by 



'l + 2tx{z + w) 



(4.34) 



(4.35) 



Ho 



(d x xf + (d xX Y 



(4.36) 



V 1 - 2tt(z + w) 

The inverse transformation of Eq. J4.34[) . which relates 
$l and to xl and XR has the same form and it is 
obtained simply by replacing k by 1/k. 

The Hamiltonian Ho in terms of the new fields becomes 

v 
2 

where the "dimensionless velocity" v is 

v = v/l -Air 2 {z + w) 2 (4.37) 

which can be absorbed in a suitable rescaling of the x 
coordinate, x — * x^/v. Notice that the parameter k plays 
the same role as the stiffness K defined above, which gov- 
erned the change of the scaling dimensions at zero doping. 
Similarly, n governs the change of the scaling dimensions 
along the z = w phase boundary of the systems at finite 
hole density. Here too, the relationship between this stiff- 
ness k and the microscopic interactions is non-universal, 
and the validity of Eqs. (|4.35p and (|4.37|) is restricted to 
the weak coupling regime in which this continuum theory 
holds. Notice that, since z > and w > 0, we will always 
have k > 1. This fact will play an important role below. 

We can now use these results to determine the scal- 
ing dimension of the perturbation H pcr t along the phase 
boundary. It is straightforward to write the perturbation 
Hpcrt in terms of the new field x and its dual x : 



H 



pert 



2(z-w) 



sin( v 87t$l ) sin(v 87r$/?,) 



(z — w) 




(4.38) 
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C. Critical behavior along the phase boundary 



1. The correlation length exponent 



From Eq. (|4.38|) we find that the operator which per- 
turbs the line of fixed points along the phase boundary at 
finite density, 7i pcrp , involves two operators whose scal- 

2 

ing dimensions are 2k > 2 and — < 2 respectively, since 

K 

k > 1. Thus this operator becomes more relevant along 
the phase boundary, away from the KT point, which in 
this language has k = 1. In fact, if we neglect the ef- 
fects of the irrelevant operator (which is a safe thing to 
do only away from the KT point since its only important 
effect is a finite renormalization of k) we see that the ef- 
fective theory in the vicinity of the phase boundary is a 
sine-Gordon theory for the dual field \. Since the scal- 
ing dimension of the relevant operator is 2//t, it follows 
that, away from the KT point, the correlation length £ 
diverges as the phase boundary is approached as 



£ ~ \z - w\ 



1 



2-1 

K 



2(«-l) 



(4.39) 



Thus, the correlation length exponent decreases (from 
infinity!) along the phase boundary away from the KT 
point. It is apparent from the form of the perturbation 
that away from the KT point there is simple scaling, 
up to contributions of strictly irrelevant operators. On 
the other hand, as k — > 1 the relevant operator becomes 
marginally relevant and the irrelevant operator becomes 
marginally irrelevant. Thus, as K —* 1 we should ex- 
pect logarithmic corrections to scaling, and a complex 
crossover near k = 1. 



2. The columnar order parameter and its susceptibility 

On the other hand, we noted above that the dimen- 
sion of the columnar order parameter operator, i(Oi o + 
0-i,o) remains fixed at the KT value of Ai o = 1/8. 
Hence for this operator we find rjifi = 1/4. On the other 
hand from scaling we know that the susceptibility expo- 
nent obeys the scaling relation 71 q = (2 — rj)v, where v 
is given by Eq. (|4.39|) . Hence 



7k 



7i,o 



8(k- 1) 



(4.40) 



is the susceptibility exponent of the columnar order pa- 
rameter, which also increases along the phase boundary, 
even though 71,0/f = 7/4 along the whole phase bound- 
ary (provided the transition remains continuous!.) 



3. The orientational order parameter and its susceptibility 

We can use the operator identifications to look at the 
behavior of the orientational order parameter which we 



saw above is the operator O2.0 of the original version of 
the theory. We also saw that this operator is a component 
of the spin-1/2 representation of su(2)2- We also found 
how it transforms. In particular, we have 



O±2.o -> Te 



/2tt „ 

tm/ — X 



(4.41) 



Along the phase boundary the scaling dimension of this 
operator is 



1 1 

A2 ' = 2^ < 2 

which it is always relevant, and 

1 

V2,0 = ~ 



(4.42) 



(4.43) 



Using once again the scaling relation 72,0 = (2— 772,0)^, WG 
find that the susceptibility exponent for the orientational 
order parameter is 



72,0 



2k- 1 
2(k- 1) 



(4.44) 



which also decreases along the phase boundary away from 
the KT point. 



D. Tricritical Point, First-Order Transition, and 
Phase Separation 

Let us now discuss how this critical line turns into a 
first-order transition at a multicritical point. In Sections 
\V\ and IV Dl we use Monte Carlo simulations to show that 
this is indeed what happens. In Section IIV Bl we used 
mean-field methods which indicated that the transition 
eventually should become first order. For this to work 
we should be able to predict the existence of a tricritical 
point along the phase boundary at which the transition 
becomes first order. It turns out that this is the case 
and that the first-order transition is triggered by an ef- 
fective attractive interaction between holes on the same 
sublattice, leading to phase separation. 

To see how this happens we need to discuss the effects 
of irrelevant operators along the phase boundary. As we 
stated above, their most important effect is a finite and 
non-universal renormalization of k away from the value 
given in Eq. (|4.35[) . We have also focused on the role of 
the operators 04,0 and Oo,i as they are both marginal at 
the KT transition. However, we also saw that one com- 
bination of these two operators remains marginal along 
the phase boundary and its coupling constant determines 
the value of k, through Eq. (|4.35|) . On the other hand, 
the other combination is the sum of a relevant opera- 
tor, cos(^/(87r/ k)x) with scaling dimension 2/k, and of 
an irrelevant operator, cos(v / 8~7tkx) with scaling dimen- 
sion 2k. The dimension of the irrelevant operator in- 
creases along the phase boundary (thus becoming more 
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irrelevant) while the dimension of the relevant operator 
decreases as k increases (thus becoming more relevant.) 

One possible scenario for a first-order transition is 
found by noting that as k — > oo, the dimension of the 
relevant operator vanishes, and the "thermal eigenvalue" 
Joi/v ~ * 2. Thus at this point, naturally provided this 
limit is accessible, the line of fixed points reaches a dis- 
continuity fixed point 6 - 1 - and the transition becomes first 
order. However, at this point the theory becomes patho- 
logical (as k, — ► oo, v — ► 0) and one may suspect that 
other physical effects, contained in irrelevant operators, 
may intervene before this happens. 

We have so far neglected other operators which arc 
even more irrelevant at the KT transition. For exam- 
ple, the operators Os,o and Oo,2, have dimension 8 at the 
KT point. Recall that the operator Oo,2 represents pairs 
of holes on the same sublattice. Both of these operators 
are present in any lattice problem (such as the interact- 
ing dimer model) and play no significant role at the KT 
transition (beyond a non-universal but otherwise trivial 
shift of the critical coupling) and for this reason they were 
(correctly) neglected. However, along the phase bound- 
ary the scaling dimensions of these operators change. Us- 
ing the OPE, it is easy to see that along the phase bound- 
ary both operators contain the operators (among others 
which are less important) O ,2 ~ cos(2-y/ (8tt/k)x) with 
scaling dimension 8/k, and Os,o ~ cos(2v8~7tkx) with 
scaling dimension 8k. Even though they are not explic- 
itly present in our starting theory, these operators will be 
generated under renormalization and close enough to the 
KT point, k > 1, they both are and remain irrelevant. 

However, although k also changes in a non-universal 
manner, the dependence of the dimensions with k does 
not, as it follows from the structure of the theory. Thus, 
provided the dependence between k and the microscopic 
couplings allow it, it may be possible to reach a point 
along the phase boundary at which k = 4. This will 
happen at a critical value of the coupling constant u and 
a critical value of the hole density p (or, equivalently at 
a critical value of the hole fugacity z (cf. Fig. [J) ). 

At this critical value of k, the scaling dimension of the 
operators Oo,2 becomes equal to 2, and together with 
the strictly marginal operator O^o, there are now two 
marginal operators at this point. Thus the system is at 
a tricritical point at this value of the parameters^ Past 
this point, Oo,2 becomes marginally relevant along the 
phase boundary. In this regime, the effective field theory 
at the phase boundary is a sine-Gordon theory for the 
field x with the marginally relevant operator Oo,2- Since 
the sine-Gordon theory in this regime is massive, it has 
a finite correlation length and since O0.2 is marginally 
relevant the correlation length along the phase bound- 
ary, which has now become a coexistence curve, has an 
essential singularity as a function of the distance to the 
tricritical point, i.e. a KT-like transition. Thus, the tran- 
sition becomes first order along the phase boundary past 
the tricritical point with a correlation length that scales 
like £ ~ e const -/ v*, where s is the distance to the tricriti- 



cal point measured along the coexistence curve. In con- 
trast, the correlation length across the phase boundary 
(below the tricritical point) exhibits conventional power- 
law scaling. Closely related scenarios for the existence 
of such tricritical points have been suggested in other 
systems, such as the extended Hubbard model in one- 
dimension^, the two-dimensional classical Ashkin- Teller 
mode l 64 ! 65 , and the dilute 4-state Potts model 6 - 6 -, which 
is a statistical system with very similar phase diagram. 

What happens as the tricritical point is reached, can 
be understood more physically by noting that at that 
point the operator Oo,2, which measures the probability 
amplitude for a pair of holes (on the same sublattice), 
becomes relevant. The relevance of Oo,2 indicates that 
holes on the same sublattice now have a strong effective 
attractive interaction, have a strong tendency to pair- 
ing and consequently phase separate. The effective field 
theory description given above corresponds to the grand- 
canonical picture, since the coupling constants arc simple 
functions of the hole fugacity. On the other hand, in the 
canonical description, i.e. at fixed hole density p, the co- 
existence curve opens up into a two-phase region: there 
is phase separation between hole-poor regions with local 
columnar dimer order and hole-rich regions. The jump 
in the hole and dimer densities (as well as in the order 
parameters) across the first-order transition is governed 
by the correlation length at the coexistence curve. Thus, 
close to the tricritical point the jump in the densities 
{i.e. the width in density of the two-phase region) has 
the scaling form Ap ~ and therefore vanishes with an 
essential singularity as the tricritical point is approached. 
Similar scaling behavior applies to the discontinuity of 
the columnar and orientational order parameters across 
the two-phase region. 

In the subsequent sections we will give further evidence 
for the nature of the phase transitions in this system, in- 
cluding the first-order transition, using Monte Carlo sim- 
ulations in the canonical and grand-canonical ensemble. 



V. MONTE CARLO SIMULATIONS 

We now employ Monte Carlo simulations to map out 
the phase diagram of the doped quantum dimer models 
at their generalized RK points. This approach is comple- 
mentary to the analytic approach of Sections IIII1 IIV Al 
and II V Bl and of Appendix [A"l We first introduce (sub- 
section lV A[) a canonical Monte Carlo algorithm for inter- 
acting dimers. In subsection IV Bl we apply this method 
to the case of the interacting fully packed classical dimer 
model on the square lattice, a system that has been stud- 
ied recently in some detail ) 14 i 16 ' 17 i 18 and study its phase 
transition. In subsections IV CI and IV Dl we consider the 
case of the doped dimer model at low doping and map 
out the critical line, verifying the theoretical scenario dis- 
cussed in Section IIV Bl In subsection IV Dl we combine 
the cluster algorithm with conventional grand-canonical 
moves that permit us to determine the first-order transi- 
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tion line as a function of dimcr fugacity z<i and interaction 
strength u. We also estimate the location of the multi- 
critical point discussed in Appendix lAl and Section llVBl 
using data from both canonical and grand-canonical sim- 
ulations. 



A. Algorithm for classical interacting dimers 

At high dimer coverage (low doping), conventional 
Monte Carlo algorithms become very inefficient. On the 
other hand, in Ref. [H, it was demonstrated that a geo- 
metric cluster algorithm (GCA) can work efficiently for 
dimers that only have a repulsive hard-core interaction. 
We briefly summarize this algorithm here. The overlap 
of two hard-core dimer configurations generates a transi- 
tion graph. This graph consists of disjoint subgraphs of 
dimers alternating between the two configurations. In the 
presence of holes, there are two possible types of graphs: 
an open graph which always terminates on a hole or a 
closed loop. Any Monte Carlo move corresponds to a 
transition graph of the initial and final configurations. 
In the geometric cluster algorithm, the two subgraphs 
are related by a global lattice symmetry. The algorithm 
obtains long transition graphs with minimal overhead: 
moves are never rejected, and each dimer encountered 
during the construction of the graph participates in the 
move. The construction proceeds as follows^ First, a 
"seed" dimer and a symmetry axis are chosen at random. 
The seed dimer is reflected with respect to the symmetry 
axis, and if it overlaps with other dimers these are re- 
flected as well. This proceeds in an iterative fashion until 
there are no more dimer overlaps or, equivalently, when 
an open or closed graph has been formed. On the square 
lattice, the algorithm is crgodic if we allow both diagonal 
and horizontal-vertical axes passing through sites of the 
lattice. The first choice allows to change the numbers of 
horizontal and vertical dimers, whereas the second one 
permits to move through the different winding number 
sectors. Transition graphs generated by the algorithm 
are symmetric with respect to the symmetry axis, and 
cross it at most twice. 

We now extend this approach to dimers with addi- 
tional interactions by exploiting the generalized geomet- 
ric cluster algorithm proposed by Liu and Luijteni 24 i 67 
Now, in a single cluster move, multiple transition graphs 
and/or open graphs are formed simultaneously while re- 
taining the rejection-free character of the algorithm. This 
is achieved by also reflecting dimers that do not overlap, 
with a probability that depends on the dimer-dimer cou- 
pling. When a dimer i, located at r° ld , is reflected to a 
new position ff cw , there are two classes of dimers that 
interact with dimer i: a) dimers which interact with it 
before it is reflected and b) dimers which interact with i 
after it is reflected. Dimers j, located at positions r*j, 
that belong to any of the two classes are included in the 
cluster (i.e., will be reflected with respect to the symme- 



try axis) with a probability 



SUij 



(5.1) 



where 5U i3 = V{\rf cw -rj\) -V{\r^ ld -r 3 \) and V(r) rep- 
resents the interaction between two dimers at a separa- 
tion r. Thus the cluster addition probability for dimer j 
depends solely on the energy difference corresponding to 
a change in relative position of i and j. In the limit 
of a pure hard-core repulsion, this generalized geometric 
cluster algorithm reduces to the original GCA. 

The GGCA applies only to simulations in the canon- 
ical ensemble. To perform Monte Carlo simulations 
in the grand-canonical ensemble, we alternate the clus- 
ter moves with conventional grand-canonical Metropolis 
moves, consisting of insertion and deletion attempts of 
single dimers. 



B. 



Zero doping: Kosterlitz-Thouless transition to a 
columnar valence-bond crystal 



According to the theoretical study of Section [IV Al and 
also from the results of Refs. [I^fl6lfl7l . we expect to find a 
Kosterlitz-Thouless transition at zero doping, as a func- 
tion of the dimer interaction. To detect and locate this 
transition, we exploit the fact that there is an ordered 
phase in the large-u region and define columnar and ori- 
cntational order parameters, 



C(r) = h(r)-n.(r + e,)] 



(5.2) 



R(r) = n x (r)n x (r + e y ) - n y (r)n y (r + e x ) , (5.3) 

where rii(r) denotes the dimer density at r. (C(r)) is 
non- vanishing only in a columnar-ordered phase, pro- 
viding a signature of translational symmetry breaking 
(with a four- fold degeneracy), whereas (R(r)) measures 
the breaking of invariance under 7r/2 rotations. In terms 
of the most relevant operators of the effective theory of 
Section TlV A| using OPE, we make the following identi- 
fications: 



C ~ -(Oi,o + 0-i,o) 
R ~ ^(O 2 ,o + O„ 2 ,o) 



(5.4) 
(5.5) 



Having a proper correspondence between the effective 
theory described in Section IIV Al and the microscopic 
order parameters (|5.3[) . we may verify our predictions. 
Since the conventional fourth-moment ratio (directly re- 
lated to the Binder cumulant££) of the order parameter 
generally does not show a well-defined crossing at a KT 
transition, we instead use a scaling function of the form of 
the spontaneous staggered polarization of the six-vertex 
models which maps on the same vertex operator as R 
in the Coulomb-gas representation, and is in the same 
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universality class. Keeping only the most relevant terms, 
this function has the form, 



(R(u)) = a 



1 



(5.6) 

From a careful nonlinear least-squares fit to the numerical 
data outside the finite-size regime (cf. Fig. [3|) we obtain 
u r = 1.508 ±0.003. 
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FIG. 3: (color online) The orientational order parameter R 
in the undoped case: it vanishes for u < u c and has an es- 
sential singularity at u c . The curves represent Monte Carlo 
data for different lattice sizes, interpolated via multiple his- 
togram reweighting. Inset: data collapse, with a least-squares 
fit to the exact scaling function for the staggered polarization 
operator of the six-vertex model (|5.6[) related to R by a uni- 
versality mapping as discussed in the text. 





FIG. 5: (color online) Fourth-order amplitude ratio of the 
orientational order parameter R in the undoped case. In con- 
trast to the amplitude ratio of the columnar order parameter 
C (Fig. |4]) , this quantity exhibits a strong hnite-size depen- 
dence, leading to an effective crossing point that can be used 
to locate the transition point. 



with M = C, R. The behavior of Qc, shown in Fig. 2J 
is similar to what is expected for the XY model, namely 
a collapse of all curves in the critical low-u phase and no 
well-defined crossing of curves for different system sizes. 
In contrast, Qr (Fig. [5]) is found to exhibit such strong 
finite-size effects in the critical phase that its behavior al- 
most resembles that of a regular continuous phase transi- 
tion. This anomalous behavior explains why the crossing 
point of the curves for different system sizes could be 
exploited to obtain an accurate estimate of the critical 
couplin g 14 ' 16 . In the figures presented in this section, 
the multiple histogram reweighting method^ has been 
used to interpolate all data obtained at different values 
for the coupling parameter. This allows us to accurately 
locate crossing points and extrema in the curves. 

Alet and coworker a 14 i 16 and Poilblanc and coworkers^ 
used transfer-matrix calculations and Monte Carlo sim- 
ulations to study the critical behavior of the undoped 
system for u > (attractive dimer interactions), whereas 
Castelnovo and coworkers^ used transfer-matrix meth- 
ods to study primarily the u < ("repulsive") regime. In 
addition, in Ref. [IH the doped interacting dimer model 
was also briefly studied for low doping by means of nu- 
merical transfer- matrix techniques. All these results are 
consistent and complementary to those presented in the 
following subsection. 



FIG. 4: (color online) Fourth-order amplitude ratio of the 
columnar order parameter C in the undoped case. The curves 
for all system sizes essentially coincide for the entire critical 
phase (u < uc), as expected for a KT transition. 

For completeness, we also investigate the behavior of 
the fourth-order amplitude ratios Qm = (M 2 ) 2 / (M) 4 



C. Low doping: Line of fixed points 

We now proceed to the low-doping regime. We per- 
form simulations in the canonical ensemble, for cou- 
plings near the critical region and for hole densities 
p h = 0.004,0.01,0.02,0.04,0.06. Figure © shows, for 
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FIG. 6: (color online) Classical susceptibilities of (a) the 
columnar order parameter C and (b) the orientational order 
parameter R, for ph = 0.06. The magnitude and position of 
the maxima for different lattice sizes L can be used to extract 
the critical exponents 7 and v. 
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FIG. 7: (color online) Fourth-order amplitude ratios of (a) the 
columnar order parameter C and (b) the orientational order 
parameter R, for ph = 0.06. The clear crossings of the curves 
for different system sizes indicate a regular continuous phase 
transition. 



Ph = 0.06, the susceptibilities of the columnar and ori- 
entational order parameters, C and R, followed by the 
corresponding fourth-order amplitude ratios in Fig. [JJ 

According to the predictions of Sect ion HV B[ we expect 
that the columnar order parameter C, which maps onto 
the O1.0 effective operator in the scaling limit, will re- 
tain its scaling dimension 1/8 along the critical line that 
emerges from the KT transition point for low hole doping 
and which constitutes the phase boundary between the 
dimer-holc liquid and the columnar solid phases. This 
constant value of the scaling dimension of C is the most 
salient signature of this critical line: The scaling dimen- 
sions of all the other operators change continuously along 
the phase boundary. 

To test this prediction, we extract the anomalous di- 
mensions rjc and r]R (which are equal to twice their scal- 
ing dimension) by means of finite-size scaling. The max- 
imum of the susceptibility scales as ^ max ~ L 2 ~ v + ■ • • . 
Subleading scaling contributions are omitted in the fit- 
ting expression since the results for sufficiently large lat- 
tice sizes satisfy simple scaling [cf. Figs. [S] and [§]. 

The correlation-length exponent v can be extracted 



from the slope of the fourth-order amplitude ratios of 
both order parameters at the critical point. Here, in- 
stead, we obtain it from the scaling behavior of the loca- 
tion of susceptibility maximum, which scales as it x max = 
u c + const. L- 1 ^ + ■■■ (cf. Figs. |8(b)| and |9(b)[). The 



critical coupling u c , in turn, is obtained from the cross- 
ing points of the fourth-order amplitude ratio (Fig. [7]). 

By repeating this procedure for different hole densities, 
we find rjc and 77^, as well as v as a function of ph- We 
note that for the undoped case the (logarithmic) finite- 
size corrections are so strong that the anomalous expo- 
nents are very difficult to determine. By including sub- 
leading corrections to the susceptibility expressions at the 
KT transition, we find strong indications that i]c = 1/4 
and r]n = 1, satisfying the established theoretical de- 
scription of this KT transition. The density dependence 
of the exponents is shown in Fig. 10(a) Clearly, they be- 
have very differently: For the columnar parameter C, its 
anomalous dimension remains unchanged and equal to 
1/4, while for R it decreases monotonically from 1 to 1/4 
where the transition is expected to become first-order, 
according to the scenario presented in Section HVB1 The 
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FIG. 8: (color online) Finite-size scaling analysis of the sus- 
ceptibility of the columnar order parameter C for ph = 0.06: 
(a) Scaling of the height of the maximum as a function of sys- 
tem size L. (b) Scaling of the position of the maximum. The 
critical value u c is obtained from the crossing of the fourth- 
order cumulants (cf. Fig. ■ The same power-law behavior 
of the susceptibility maximum is found for all hole densities 
p h < 0.06. 
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FIG. 9: (color online) Finite-size scaling analysis of the sus- 
ceptibility of the orientational order parameter R for pn = 
0.06: (a) Scaling of the height of the maximum as a function 
of system size L. (b) Scaling of the position of the maxi- 
mum. The critical value u c is obtained from the crossing of 
the fourth-order cumulants (cf. Fig. [TJl . Unlike the results for 
the columnar order parameter (Fig. [8} the power-law behav- 
ior of the susceptibility maximum is is now dependent on the 
hole density, consistent with an anomalous scaling exponent 
that varies along the phase boundary. 



results from our Monte Carlo simulations are thus con- 
sistent with the predictions we made in our theoretical 
analysis. 

The evolution of the correlation length exponent along 
the phase boundary is shown in Fig. 10(b) This expo- 



nent behaves qualitatively as predicted for finite doping, 
i.e., it exhibits a monotonically decreasing behavior along 
the line of fixed points. A direct quantitative comparison 
to the field-theoretical prediction is not possible, since 
the simulations are performed in the canonical ensem- 
ble. Even though the dimer density could be mapped to 
a dimer fugacity, the field theory assumes a fixed hole 
fugacity. In addition, the evolution of the exponent in 
Fig. 1 10 (b) | is slower than predicted because in the simu- 
lations we do not approach the phase boundary perpen- 
dicularly in the simulations, which leads to an effective 
exponent v that is smaller than the one computed from 



the scaling dimension of the relevant operator in the field 
theory. 

In Fig. [11] we summarize our results for the location of 
the phase boundary between the dimer-hole liquid phase 
and the columnar solid phase for hole densities ph < 0.06. 
The behavior of the critical line for ph — > is consistent 
with its expected scaling behavior ph oc £,~ 2 (ph = 0) 
which is based on simple dimensional analysis, with the 
only length scale of the problem being the correlation 
length of the undoped problem. 



D. High doping: 



First Order Transition and Phase 
Separation 



Beyond the multicritical point predicted in Sections 
IIIII and IIVBI (see also Appendix [XJ , we expect a first- 
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FIG. 10: (a) Evolution of the anomalous dimensions 77c (lower 
curve) and r\n (upper curve) for the columnar and orienta- 
tional order parameters, respectively, along the critical line 
for hole densities 0.004 < ph < 0.06. The transition will be- 
come first-order when the two curves cross, (b) The inverse 
correlation- length exponent \jv as a function of doping along 
the critical line for hole densities 0.004 < ph < 0.06. Up to 
small systematic deviations discussed in the text, the observed 
trend agrees with the scaling predictions. 



order transition line in the zj-u phase diagram, where 
Zd denotes the dimer fugacity. This line separates the 
crystalline from the liquid phase. According to our sce- 
nario, an entropic attraction between holes on the same 
sublattice becomes marginal and leads naturally to phase 
separation between a hole-rich liquid phase and a hole- 
poor columnar crystalline phase. Since the multicritical 
point is characterized by a marginally relevant opera- 
tor, complicated crossover will be observed in the first- 
order transition region in the vicinity of this points In 
addition, close to the multicritical point, the first-order 
transition will be very weak, with a discontinuity that 
vanishes with an essential singularity as a function of the 
distance to the multicritical point along the phase coex- 
istence curve. At this transition, all observables, such as 
the latent heat, should vanish in a similar way, making 
the numerical study of the transition close to the multi- 
critical point particularly difficult. 

To confirm the existence of the discontinuous transi- 
tion we perform grand-canonical Monte Carlo simula- 
tions with singlc-dimcr insertions and deletions (alter- 
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FIG. 11: Phase diagram at low doping, with the critical line of 
fixed points separating the dimer-hole liquid phase [u < u c (x)\ 
from the columnar solid [u > u c (x)]. 



nated with canonical geometrical cluster moves to accel- 
erate the relaxation of the configurations), for couplings 
u = 3.0, 3.5, 4.0, 5.0, 6.0, as a function of the dimer fugac- 
ity Zd- Figure 12(a) shows the dimer density as a function 
of dimer chemical potential Hd- Although with increasing 
system size a jump in the dimer density develops, it does 
not become very pronounced. However, consideration of 
the heat capacity Cy [Fig. 12(b)| confirms the presence of 
a single phase transition at fixed coupling constant u, as 
Cy exhibits a peak at a chemical potential that matches 
the location of the jump in pd- We emphasize that this 
classical heat capacity is not the heat capacity of the 
(2 + l)-dimensional QDM. Indeed, Cy does not have 
any physical meaning in terms of the ground-state wave 
function that we are considering, because the ground- 
state energy cannot be changed through variation of the 
parameters u or pid- Another confirmation of the phase 
transition is obtained from the susceptibilities of the ori- 
entational [Fig. 13(a)| and columnar [Fig. 14(a)| order 
parameters, \R an d Xc, respectively. Both quantities 
exhibit a peak at a chemical potential that approaches, 
with increasing system size, the location of the peak ob- 
served in Cy. 

To confirm the nature of the phase transition, we con- 
sider the scaling of the peaks in Cy, XRi an d Xc- For a 
first-order transition, these quantities should exhibit a 5- 
function singularity in the thermodynamic limit or equiv- 
alently, for finite systems their peaks must scale with the 
lattice size in a finite system. We find that the heat- 
capacity peak, apart from a constant background, indeed 
scales with the lattice size L 2 for the range of system sizes 
(up to L = 140) that we considered. This is supported 
by the behavior of the system-size dependent maxima in 
Xr and xc, see Figs. [l~3l and [l~4l which both scale as L 2 , 
indicating that r) = 0. In addition, all local order param- 
eters must develop a discontinuity at the transition point 
as the system size increases. Whereas the jump in the 
dimer density [Fig. 12(a)| is not very sharp, the jump 
in the columnar and oricntational order parameters, C 
and R, is quite pronounced already for the system sizes 
studied here, see Figs. 15(a) and |15(b)| Furthermore, 
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FIG. 12: (color online) (a) The pd-fJ-d equation of state for 
coupling u = 3.5 and different lattice sizes. As the system 
size increases, the jump in the dimer density at the first-order 
phase transition gradually becomes more pronounced, but re- 
mains relatively weak, (b) The heat capacity of the classical 
interacting dimer model for the same coupling strength. We 
find that the peak scales as C max ~ Co + L 2 , providing strong 
evidence for a first-order transition at u — 3.5. 



the location the order-parameter jump provides a good 
indication of the transition point. 

The strongest evidence, however, is provided by the 
fourth-order amplitude ratio of the density, Q p . At a 
first-order transition, this quantity displays a specific be- 
havior, as discussed in Ref. [toI . In particular, the po- 
sitions of two minima observed in Fig. 16(a) approach, 
in the thermodynamic limit, the densities of the two co- 
existing phases. Outside the coexistence region, Q p ap- 
proaches a limiting value 1/3, characteristic of Gaussian 
fluctuations. This type of behavior is not found at a con- 
tinuous phase transition, and should be considered as a 
strong indicator for the occurrence of a first-order transi- 
tion. This is particularly important since the very weak 
nature of the first-order transition makes it impossible 
to unambiguously confirm the existence of a double peak 
in the histograms of the internal energy for the system 
sizes that wc considered. Whereas the first-order transi- 
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FIG. 13: (color online) (a) Susceptibility of the orientational 
order parameter R for coupling u — 3.5 and different lat- 
tice sizes, (b) Asymptotically, the maxima of \R scale with 
the lattice size, confirming the first-order nature of the phase 
transition. 



tion becomes more pronounced at higher couplings, and 
it thus should become easier to distinguish the two peaks 
in the energy histogram, in practice those simulations are 
seriously hampered by the very large relaxation times en- 
countered for large dimer interactions. 

By repeating the analysis presented here for different 
couplings, and estimating the coexistence chemical po- 
tential from the convergence point observed at the order- 
parameter discontinuity (cf. Fig. ll5|) . wc derive the phase 
diagram in the zj-u plane, see Fig. |16(b)| 



VI. EFFECTS OF REPULSIVE HOLE-HOLE 
INTERACTIONS NEAR THE FIRST-ORDER 
TRANSITION REGION 

We now discuss briefly the effects of additional inter- 
actions near the coexistence curve. It is clear that ad- 
ditional interactions near the first-order transition line 
should stabilize more complex ordered inhomogeneous 
phases. The simplest interaction that competes with the 
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FIG. 14: (color online) (a) Susceptibility of the columnar or- 
der parameter C for coupling u = 3.5. (b) The maxima of 
\c as a function of the linear lattice dimension. As for \u 
(Fig. |13(b)[ ), the maxima scale with lattice size, supporting 
the occurrence of a 5-function singularity in \c in the thermo- 
dynamic limit, as expected for a first-order phase transition. 



tendency of holes to pair and phase separate from the 
crystal is a weak nearest-neighbor hole-hole repulsion Vh- 
The addition of such an interaction to the classical dimer 
model would lead to an additional energy cost for homo- 
geneous and isotropic clusters of holes. Remarkably, for a 
range of dimer interactions u, this energy cost leads to the 
formation of commensurate hole stripes with period 3, in 
a region in the phase diagram located between the dimer- 
columnar crystal and the hole-dimer liquid. For general 
values of u, Vh and hole densities one expects a com- 
plex phase diagram, most likely similar to what is found 
in theories of commensurate-incommensurate transitions, 
which we do not exp lore here in detail but are discussed 
in Refs. [llz21zlzl 



This phase can be thought of as the ground-state wave 
function of a quantum Hamiltonian constructed using the 
approach described in Section[H] The quantum Hamilto- 
nian that leads to the prescribed wave function includes 
a generalized form of the hole-related Hamiltonian of 



FIG. 15: (color online) (a) The columnar order parameter C 
as a function of the dimer chemical potential (id for coupling 
u = 3.5. The rapid increase near the transition becomes more 
pronounced as the thermodynamic limit is approached. The 
dashed line shows the approximate position of the transition 
point, (b) The orientational order parameter R as a function 
of the chemical potential fj,d for coupling u = 3.5. It also 
shows a rapid increase similar to the behavior of C, at the 
same chemical potential. 
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H — Hd + thoic 2_j 
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R c h' 



^)(^|+/-rV|cr)(cf 



(6.1) 



where R C h and R c h> denote the number of pairs of holes 

formed in the corresponding configurations and 
(cf. Fig. [T7]). More specifically, the ground-state wave 
function is 



\GnI) = 



1 



(6.2) 
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FIG. 17: (color online) A particular hopping process real- 
ized in the Hamiltonian Eq. (|6.1|l . The potential terms that 
are present in the Hamiltonian depend on the number of addi- 
tional pairs of holes that are formed after the hopping process. 
In the process shown above, this number is 1. 
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FIG. 16: (color online) (a) Fourth-order amplitude ratio for 
the dimer density, at coupling it = 3.5. The exhibited behav- 
ior is typical of a first order transition.—, with two minima 
that in the thermodynamic limit approach the coexisting den- 
sities, (b) Phase boundary between the dimer liquid and the 
columnar solid. For low couplings, the transition is continu- 
ous (solid line). Our best estimate for the multicritical point 
(u c ~ 2.6) is represented by the large circular dot. Beyond the 
multicritical point, the phase boundary (dashed line) is esti- 
mated from grand-canonical MC simulations at fixed dimer 
fugacity Zd- 



This wave function has a counterpart in the grand- 
canonical ensemble, for which the exactly solvable quan- 
tum Hamiltonian is a generalization, in exactly the same 
way as above, of Eq. (|2.8[) . 

By performing grand-canonical simulations for weak 
hole repulsions Vh < jqU in the regime of strong dimer 
attractions, u > 4.0, where the first-order transition is 
more pronounced, a hole-stripe phase was observed. In 
particular, for the couplings u = 5.0, Vh = 0.5 and a 
dimer chemical potential —1.1 < /Jd < —0.8 (i.e., be- 
tween the liquid phase /i^ < —1.2 and the columnar phase 
Md — 0.8), the hole density structure factor shows non- 
trivial peaks at (k x ,k y ) = (±27r/3, 0), (0, ±27r/3), which 



FIG. 18: (color online) The hole density structure factor S(k) 
for dimer coupling u — 5.0, hole repulsion Vh = 0.5 and dimer 
chemical potential [id = —1.0, for linear system size L = 128. 
The four peaks at (±27r/3, 0) and (0, ±27r/3) correspond to 
the ordered configuration of hole-stripes shown in Fig. 1191 



sharpen as the lattice size increases (see Fig. [T5|) . A snap- 
shot of the ordered phase (Fig. [19)1 illustrates how the 
holes order in stripes with a period of three lattice spac- 
ings, whereas the dimers are still ordered in a columnar 
pattern. In this way, the holes minimize the effect of the 
hole-hole repulsions and the dimers simultaneously maxi- 
mize the effect of the attractive dimer-dimer interactions. 
The same pattern is also found for u = 4.0,4.5,6.0, in 
similar regimes for the hole-hole interaction. 

More generally, we expect that, as the liquid phase is 
approached in the regime of strong dimer couplings, a 
sequence of hole-commensurate phases will be stabilized, 
leading ultimately to incommensurate phases next to the 
liquid phase. The formation of this phase diagram is 
similar in spirit to the ones discussed in Refs. I7l1l72ll73ll74l . 
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FIG. 19: Snapshot of a part of a typical ordered configuration 
that appears at the couplings u = 5.0, Vh = 0.5, fid = —1.0 
and linear system size L — 128. Holes prefer to form com- 
mensurate stripes of period 3, so as to minimize the effect of 
the weak hole repulsions. 



VII. ELEMENTARY QUANTUM 
EXCITATIONS OF DOPED QUANTUM DIMER 
MODELS 

In the preceding two sections we discussed the proper- 
ties of the ground-state wave functions and the behavior 
of equal-time correlation functions of several operators of 
physical interest. There we used extensively the connec- 
tion that exists for these type of wave functions between 
the computation of equal-time correlators of local opera- 
tors and computations of similar objects in the equivalent 
two-dimensional classical statistical mechanical system of 
interacting dimers and holes. In this section we will be 
interested in the spectrum of low lying excitations which 
is inherently a quantum mechanical property. 

Unfortunately, as it usually the case in QDMs,— all we 
know is the ground-state wave function. The low-lying 
excitation spectrum is not known exactly but it can be 
computed approximately using the variational principle. 
This is the single mode approximation (SMA), which is a 
useful tool for studying the excited states of many body 
systemSi 75 ' 76 i 77 ' 78 i 79 It is particularly useful in the case of 
quantum dimer models at their RK points due to the fact 
that the exact ground-state wave function is known ex- 
actly. The computation of the low lying collective modes 
in the QDM was done by Rokhsar and Kivelson.— Al- 
ternatively, one can describe qualitatively the low-lying 
spectrum using the effective field theory of the quantum 
dimer models (and their generalizations) at criticality, 
the quantum Lifshitz model of Ref. [Til 

In this section we will consider only the SMA spectrum 
in the dimer-hole liquid phase. Similar calculations can 
be done in the phase with long-range columnar order. 
In the dimer-hole liquid phase the equal-time correlation 
function of the hole density operator, i.e. the one-body 



density matrix, approaches a constant at long distances. 
Thus the wave function for this phase exhibits a Bose 
condensate of holes. Since the holes are charged, this is 
a charge Bose condensate. To determine if it is a super- 
fluid (or more precisely a superconductor) it is necessary 
to show that it has a finite supcrfiuid stiffness, i.e. a crit- 
ical velocity. This can be determined form the spectrum 
of density fluctuations and hence from the spectrum of 
collective modes. It will turn out that, in spite of the 
more correlated nature of the wave functions we consider 
here, the result will be similar to that of Rokhsar and 



Kivelson? 



no superfluid stiffness. Given the more 



general structure (while still local) of the wave functions 
we study here, we conclude that wave functions associ- 
ated with Hamiltonians satisfying the RK condition in 
general do not describe superfluid states. 

We begin by summarizing the SMA procedure, focus- 
ing on doped quantum dimer models. Firstly, one must 
know the exact ground state of the system |0) and the 
type of excitations which saturate the frequency sum 
rule. In our case, there are two candidates: the dimer 
density and the hole density excitations. Since it follows 
from a variational principle, the SMA provides a proof of 
existence only for gaplcss excitations but not for gapped 
ones. The energy of an excitation created by an operator 
with wave vector k, which we will denote by /5k, acting on 
the ground state is bounded from above as follow o 75 ' 76 ! 77 



E < 



/(k) (0|[p(-k),[W,p(k)]]|0) 



*(k) 



(0|p(-k)p(k)|0) 



(7.1) 



where /(k) is the "oscillator strength" and s(k) is the 
structure factor (i.e. the equal-time correlation function) 
for the operator p(k). 

In the case of doped QDMs at their RK point we 
know their ground-state wave functions exactly and they 
have (by construction) zero ground-state energy, Eq = 0. 
Thus the excitation spectrum must be positive. The sys- 
tem will have gapless excitations if E^ — Eq vanishes at 
some wave vector. It is worth noting that there are two 
distinct ways for the SMA bound to vanish close to some 
k = ko. One way is if /(k) vanishes at krj. This can 
happen only if the commutator [H, /5(k,)] vanishes. This 
means that /o(krj) is a conserved quantity. The other way 
occurs when s(k) becomes infinite at ko. This is a signa- 
ture of a nearby density-ordered state like the columnar 
dimer crystal we found in the phase diagram of the dimer 
models under study, shown in Fig. [5] 

In the following, we will use the following operator def- 
initions. For dimers, <x| (r) denotes the dimer density op- 
erator and it takes the values ±1 if a dimer is present 
or absent at the link which begins at the position r and 
has direction a = x,y. For holes, <r h (r) denotes the hole 
density operator and it takes the values ±1 if a hole is 
present or absent at the position r. 

Details of the derivations of the SMA oscillator 
strength functions for both models are given in Appendix 
|B1 Here we just quote the main results. 
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A. The fixed hole density model 

1. Hole density excitations 

For the fixed hole density model the SMA oscillator 
strength function /(k) for hole density excitations is 
given by 

/(k) = 4t Uolc q 2 (7.2) 

for k = (tt, 7r) + q. 

From the results of the Section [TTT] and Appendix El 
we may conclude that: For k = (tt,tt) + q, the struc- 
ture factor near (tt, tt) scales like s^'^k) cx q i^--i and 

thus, s^'^O) is a constant. Thus E^' n) - E cx q 2 . 
For k = (0, 0) + q, the results from the following section, 
which formally hold only for low density of dimers, lead 
to the conclusion that s(°'°)(k) cx and s(°'°)(0) 

is again a constant but has strong oscillatory behavior 
in real space. Thus, from the correlation function we 
computed in and Appendix [A] we conclude that these 
excitations are also quadratic in the momentum q, i.e., 
Ek — Eq < q 2 . This result is consistent with the com- 
pressibility argument, given in Ref . HI which would also 
give quadratic dispersion in this case (keeping in mind 
that in this case the compressibility is infinite when the 
system is doped (constant number of holes)). However, it 
is important to stress that the compressibility argument 
is a stronger condition because our calculation of the cor- 
relation function is legitimate for low dimer densities. 

2. Dimer density excitations 

For the fixed hole density model the SMA oscillator 
strength function /(k) for dimer density excitations is 
given by 

/(k) = /dimer-fiip(k) + /holc(l)(k) (7.3) 

Close to the wave vector Qo = (tt, tt) with k = Qo + q, 
both oscillator strengths /dimcr-flip(k) and fhoie(i) vanish 
quadratically (cf. Appendix B) and more specifically, 

/dimer-flip(q) = 8t<? 2 (7.4) 

/hole(i)(q) = -^hoic<7 2 (7.5) 

Given the fact that the dimer density structure factor 
is a constant at Qo = (tt, tt) and combining Eqs. (|B12j) 
and (|7.3[) . we may conclude that there are gapless dimer 
density excitations at Qo = (tt,tt). 

In addition to these results, we may have additional 
branches of dimer density excitations, specially when 
there is a divergence of the structure factor at some wave 
vector Qo which would would be an indication of dimer 
order at a nearby phase. In particular, this happens at 
the phase boundary between the hole-dimer liquid phase 
and the columnar-ordered crystalline phase. 



B. The fixed hole fugacity model 

As above, we can study the hole density excitations of 
the second model, the grand-canonical model, in which 
the hole density is not fixed, but is determined by a pa- 
rameter z which plays the role of a hole fugacity in the 
wave function. By the fact that there is no conservation 
of the number of holes is explicitly broken, we can predict 
that the numerator /(k) will vanish only at Qo = (tt, tt). 
This happens because of the bipartite lattice symmetry 
that enforces the number holes on each sublattice to be 
equal. For the fixed hole fugacity model the SMA oscil- 
lator strength function /(k) for hole density excitations 
is 

/hoic(2)(k) = 4i 

pairing 

(2 + cos(kj;) + cos(k y )) (7.6) 

The formula shows that the numerator /hoie(2) (k) van- 
ishes only at the wave vector Qo = (tt, tt) quadratically, 
as expected. The dimer density excitations are not gap- 
less (in the sense that the numerator does not vanish for 
any Qo) because the dimer- flip contribution vanishes at 
(tt, 7t), (0,vt) and the dimer-breaking term gives a con- 
stant contribution. This is expected due to the viola- 
tion of the dimer number conservation. It is important 
to stress that the behavior of the structure factor is ex- 
actly the same as before (in the fixed hole-density model). 
The reason is the equivalence of the canonical and grand- 
canonical ensemble in the thermodynamic limit for clas- 
sical systems. 

It is certainly worth noting that, although the two 
ground-state wave functions of Eqs. (|2.7[) and (|2.9[) have 
the same ground-state physics, due to the essential equiv- 
alence of the classical canonical and grand-canonical en- 
sembles in the thermodynamic limit, the nature of their 
quantum elementary excitations is drastically different. 

VIII. CONCLUSIONS 

In this work we have constructed generalizations of the 
quantum dimer model to include the effects of dimer cor- 
relations as well as finite hole doping in the wave function. 
Throughout we considered the case of bosonic (charged) 
holes and neglected the fermionic spinons. We have 
constructed generalized RK Hamiltonians whose ground- 
state wave functions describe the effects of (attractive) 
dimer correlations and finite hole doping. We have dis- 
cussed the rich phase diagram and critical behavior of 
three doped interacting quantum dimer models at their 
RK point using both analytic methods and numerical 
simulations. We have shown that the ground-state wave 
function embodies a complex phase diagram which con- 
sists of dimer-hole liquid and columnar phases separated 
by a critical line with varying exponents, ending at a 
multicritical point with a Kosterlitz-Thouless structure 
where the transition becomes first order. The critical be- 
havior along the low doping section of the phase bound- 
ary was investigated in detail and the predictions of our 
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scaling analysis were confirmed with large-scale Monte 
Carlo simulations. Monte Carlo simulations were also 
used to show that the transition between the dimer-hole 
liquid and the columnar solid does indeed become first 
order, and to estimate the location of the multicritical 
point and of the first-order phase boundary. 

In the high-doping regime, near the first-order transi- 
tion line, additional repulsive interactions among holes 
were shown to generate, at the expense of the two-phase 
region, an even richer phase diagram with phases in 
which the dimer-hole system becomes inhomogeneous. 
In the regime of strong dimer coupling, u = 5 = 0., 
and weak hole interactions, Vh = 0.5, we found a stripe 
phase with wave vectors (27r/3,0) and (0, 2ir/3). In gen- 
eral, we expect the two-phase region to be replaced by 
a complex phase diagram with a large number of com- 
mensurate and incommensurate phases. This physics 
is well known in the context of two-dimensional clas- 
sical statistical mechanics of systems with competing 
interactions. 71 ' 72 ' 73 However, it is interesting to see how 
it arises at the level of the exact ground-state wave func- 
tion of models of strongly correlated systems, particu- 
larly given the current interest on this type of phenom- 
ena in high-temperature superconductors and related 
systems^o^^^^i 

In this paper we have also presented an analysis of the 
low-lying excitations of the quantum system and found 
that, much as in the case of the Rokhsar-Kivelson quan- 
tum dimer model, the doped system is a Bose-Einstein 
condensate but not a superfluid, since the superfluid stiff- 
ness vanishes even though both dimers and holes interact. 
It is apparent that in order to render the Bose-Einstein 
condensate a true superfluid it is necessary to violate the 
RK condition which forces the wave function to have a 
local structure. The effects of violations to the RK condi- 
tion are poorly understood, and we have not investigated 
this important problem which is essential to determine 
the generic phase diagram of these models. 
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APPENDIX A: MEAN-FIELD THEORY FOR 
DIMERS AND HOLES 

In this Appendix we give details of the mean-field 
theory summarized in Section IIIII We will use Grass- 
mann variable methods to write down the partition func- 
tions for classical interacting and doped dimers. We will 
use these methods to derive a simple mean-field theory 
for this system. While mean-field theory, as it is well 
known, fails to give the correct critical behavior in two- 
dimensional systems, it offers a good qualitative descrip- 
tion of the phases and, surprisingly, even of the gross 
features of the phase diagram. In subsequent sections 
we will use more sophisticated analytic and numerical 
methods to study the phase transitions. 

1. Non-Interacting dimers at finite hole density 

The classical dimer-hole partition function can be for- 
mulated in terms of a Grassmann functional integral, ac- 
cording to the prescription of Ref . [32|. The classical parti- 
tion function of the dimer problem on any lattice which is 
defined by the connectivity matrix M and with fugacity 
z for the dimers and 1 for the holes is: 

Z dimcr = J V-qVtf ™\ + § E« Ma mvhi v] 

(Al) 

For the square lattice, My is 1 if i,j are nearest neigh- 
bor sites and zero otherwise. For the triangular lat- 
tice, Mij is 1 if i,j are nearest neighbors and also next- 
nearest neighbor sites, but only along one diagonal, i.e., 
if j = i + x + y or j = i — x — y. Also, the fugac- 
ity z ranges from to oo. When z — > 0, the system is 
filled with holes and there are very few dimers and when 
z — > oo the system approaches the close-packed limit with 
dimers. Both limits are worth of study due to the exis- 
tence of an important theorem by Hcilmann and Lieb 85 
which proves under rather general assumptions the ab- 
sence of any phase transitions with doping in this model. 
Thus, the identification of the phase in the few dimers 
limit is enough to conclude about the phase that the sys- 
tem enters when doped. In this section we will study 
this problem in the limit in which the dimers are dilute. 
In this regime a simple mean-field theory of the Hartree 
type is expected to be accurate^ Such a crude approx- 
imation should break near criticality, i.e. near the close 
packing limit. 

To proceed, we apply a Hubbard-Stratonovich trans- 
formation to the above partition function: 

(A2) 
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where J\f is an irrelevant normalization constant. We may 
also add sources for the hole density operators Jrjrf . In 
this way, the classical dimer partition function may be 
written as follows: 



(A3) 



where he have dropped the normalization constant J\f, as 
we will do in what follows. (Note that what makes the 
above problem unsolvable is the term '1' in the exponent!) 
Upon integrating out the Grassmann variables we find: 



(A4) 

In the limit z — ► we may formulate a legitimate and 
well-defined mean-field theory. To this end we rewrite 
the partition function in the following way: 

= J V<j) e -i sw (A5) 

where 

s(4>) = ^X>*- JtWyHh- Jj)-z^2H4H + i) 

ij i 

(A6) 

is the effective action. As z — > 0, we have a theory 
which has a well-defined saddle point and the perturba- 
tion around this point will be in powers of z which plays 
the role of an effective coupling constant. In this way, we 
have a very fast convergent expansion. 
The saddle-point is defined as follows: 



5S 



= 



<t>i=4> 



and we take the following equation for 



rh: 4- 



1 



(A7) 



(A8) 



In this approximation, the density of holes is given by: 
f x d In ZfH mer 



Pi = (ViVi) 



1 dS 

z dJ t 
1 

H + 1 



(A9) 
(A10) 
(All) 



So, the source Jj in terms of the density of the holes pi, 
is given by: 



Ji 



y My 



— -z^MijPj - 1 



(A12) 
(A13) 



The Legendre transform of the effective action S is: 
T(jh) = - z S{4> j {p i ),J j { Pi ))+Y J Ji{Pj)pi (A14) 

i 

ij i i 

(A15) 

We specialize now to the case of uniform hole density 
Pi = p and also for the case of the square lattice where 
the number of nearest neighbors is 2D = 4 and assuming 
that the number of sites on the lattice is N. Then, 



t(p) 



'-(ANp 2 ) + Nln(p) +N(1 - p) (A16) 



At this level of approximation ("Hartree") the equation 
of state becomes 



J = -4zp+ - -1 
P 



(A17) 



So, when the sources are set to zero, the density in terms 
of the fugacity, in the limit z — > 0, is: 



P = 



l + VT+Wz 



(A18) 



As a check of the approximation, we may expand the 
result in the region of z — > 0, where the result should be 
p ~ 1 for the hole density: 



p = 1 -Az + 0{z 2 ) 
Or, equivalently, for p — ► 1~, 



(A19) 



g =^{^- 1 ) =\{l-p) + 0{{l-pf) (A20) 

The hole density-density correlation function can be ob- 
tained in the following way by using the Legendre trans- 
form: 



(ViViVjVj) ~ (ViVDiVoV 

d 2 In Z dimer dpi 



,t\ 



dJidJj 



dJi 



and also, 



Gij = 



' d 2 T ' 

dpidpj 



(A21) 



(A22) 
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By using (|A15|) and (|A22j) . we have: 
d 2 T 



dpidpj 



pi 



(A23) 



Since Mij is a function of the distance between sites and 
vanishes except for nearest neighbors, its Fourier trans- 
form is: 

M(q) = a 2 ^e- tq ' (r '- r ^M(r i - r 3 ) 



2a 2 cos<7 Q a 

a=l 



(A24) 



Also, we set a = 1 and finally the hole density connected 
correlation function is: 



5(q) = 



1 



j? - 2zJ2 a=1 cos q a 



(A25) 



For momenta near Q = (it, tt), q = Q + p, with p small, 
it becomes 



Q(p + Q) 

where £ is the correlation length 



£- 2 + p 2 



4/; 



(A26) 



(A27) 



Finally, the connected hole density correlation function 
(the structure factor) in real space for p — > 1 is: 



G(r) 



(_l)r. 



(A28) 



If we restore the units, the correlation length is £ 



/ l-pa 2 



Surprisingly, given how crude this approximation is, 
the result of Eq. (|A28[) is consistent with the numerical 
results provided by Krauth and Moessner— for much of 



the dimer density range they studied. Significant devi- 
ations are seen only upon approach of the close packing 
regime where the classical dimer model is critical and this 
mean-field calculation fails, e.g. the correlation length di- 
verges as p — ► with exponent 1/2, given by Eq. (|A27I) . 
(the mean-field value) . The correct value of the exponent 
can be deduced from Table U and it is 1/(2 — 1/4) = 4/7 
(1/4 being the scaling dimension of the hole operator for 
the non-interacting case.) As we show in Section IIV Al 
(and in Tabled, the dimension of the hole operator grows 
from the value 1/4 for free dimers to a value of 2 at the 
(Kosterlitz-Thouless) transition to the columnar state, 
where it should exhibit an essential singularity due to 
the marginality of the hole operator. 

2. Adding interactions between dimers 

Clearly, the mean-field method for the non-interacting 
dimer-hole system can be easily extended to dimer-hole 
systems with local interactions. In the case of attrac- 
tive interactions between parallel dimers, the partition 
function of the monomcr-dimcr system should be 



d — int. 



J Vr{Dtf cxp ( ]T r ilV l + | Y M, 

V i ij 



V 



Y M n h Vi vl Vj v} Vk vl m vl 



ijkl 



(A29) 



where V = z 2 (e u — 1) for an attractive interaction of 
strength u > between parallel neighboring dimers. Mjj 
represents the coordination array of the lattice and it 
takes the value 1 when i is nearest neighbor to j and 
otherwise is zero. In the same respect, Mijki takes the 
value 1 only when i,j,k,l are arranged on a square pla- 
quette and is zero otherwise. The sums run along all 
possible lattice sites for each index. We may perform 
two Hubbard-Stratonovich transformations by introduc- 
ing the fields \i 4> an d then we have (again, dropping all 
normalization constants) : 



J 



-'d— int 



= / 2??/Dr^2?xexp 



Y Vl] - y Y 1 Xki + Y ViViVjVj (Xij + ~Mi. 



ijkl 



T>rfDifT>xD<$> exp 



- Y XnM l3kl )- x X ki ~lY WW + \ M nY X <$>3 + Y + !) 



ijkl 



= / T>\D(j)exp 



- T7 Z (Mm) x Xki -\Y fobcv + 7; M ij) l oj + y + 1] 



ijkl 



(A30) 
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Just as we did above, we introduce a set of auxiliary sources J*- defined on links ij, and J( defined on sites i. The 



Z d _ int [J*,J*] = We ^ 



partition function now reads 



where the action now is: 



(A31) 



/ / A7 



We define the densities conjugate to the fields fc, \ij as: 
d In Z d _ int 1 



E(Xy + «M (i )-^ 



5 111 2sd—int 2 ^ — -s ~ , _-| 

m y = Klx = - 2J< M w) 



<9/ x 



r 



(A32) 



(A33) 



hi 



Solving for the fields and replacing in Eq. (|A32[ ). we have finally for the effective potential or equivalcntly the Gibbs 
free energy : 

T(z, V) = — mijMijkimkt + E Ui ( \ E ^vM m ki + ) 



ijkl 



hi 



2 L j \ fc/ 



+ -My Wlj + f 



(A34) 



The ordered phase of dimers is expected to be a columnar one. So, for my = (— f ) Xi i5i J _sm + too and rij = n, the 
free energy is: 

T<kZ ^ = Vm 2 + 2Vml + 4Vm n 2 + 2zn 2 



(A35) 



In l + 8n V(— +m ) + 



In i + 8n V(— - +rn ) + 



Now, we proceed by solving two of the three extremal (saddle-point) equations: 



ST , «jr 

— = and = 

on otoq 



(A36) 



which take the explicit form 



8VmQii + Azn — 



Wm + AVn 2 



4 (V(f+mo) + §) 4(V(-f+mo) + 



l + 8n V(f +TO ) + f l + 8n( V(-f + mo) + 



4nF 



l + 8n F(-f +m ) + f l + 8n V(f +m ) + 



= 



= 



(A37) 
(A38) 



For z = V = 0, the solution of the saddle-point equa- Eqns. (|A37[) and (|A38[) recursively, expanding also in the 
tions is trivial toq = n = 1. For small z, we can solve 
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small order parameter to (this is correct close to a con- sition): 
tinuous phase transition or to a weakly first order tran- 



m 



1 - 4z - SV 



1 + 8V + 4z 



16V 2 m 2 



256V 4 m 4 



4096U 6 to 6 



A(2V + z) y {l + 8V + Az) 2 (l + 8y + 4z) 4 (l + 8V + 4zY 



+ 



2048V 



5^4 



+ 



4V y (1 + 8V + 4z) 3 (l + 8V + 4z) 5 (l + W + 4z) 



Replacing Eqs. (|A39[) . (|A40p in the free energy (|A35[) . and also expanding in the small parameter z, we have: 



where V r = tt 



r(z,iQ 

2V 



1 and 



C (z, + C 2 (z, V r )m 2 + C 4 (z, V r )m 4 + C 6 (z, V r )m 6 + ■■■ 



(A39) 
(A40) 

(A41) 
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CoO, V r ) = -2z + (8 - 2y r )z 2 + — (-5 + W r )z 3 + 0(z 4 ) 



(A42) 

C 2 {z,V r ) = V r z 2 (1 + 8V r z 2 -128V r z 3 -256{-7 + V r )V r z 4 ) + 0(z 7 ) (A43) 
C 4 (z, V r ) = -32V 4 z 7 (l - 2(19 + V r )z + 4(184 + V r {4 + V r ))z 2 + 8(1336 + V r (-232 + V r (4 + V r )))z 3 ) + 0{z 14 ) 

(A44) 

C 6 (z, V r ) = 128V; 6 z 10 (l - 4(K + 12)2 + 4(340 + V r (28 + W r ))z 2 ) + 0(z 13 ) (A45) 

I 



In the limit of very low dimcr density, z — > 0, from Eqs. 
(|A42|) - (|A45|) . there is a clear first-order phase transition 
from an empty lattice (n = 1, m = 0) to a columnar 
dimer crystal (to ^ 0) because C 2 > 0, C4 < and C*6 > 
0. In particular, in this limit, the first-order transition 
happens when: 



C 2 



Z 2 V r 



cl 

(32V 4 z 7 ) 2 
h!2V?z la 



(A46) 
(A47) 

2z 2 e u = 1 (A48) 



The condition (|A"4"8)) can be derived through a very sim- 
ple argument: In the limit of very low dimer densities, 
the non-local effects which are related to the hard-core 
dimer constraint are negligible. If we consider just a sin- 
gle plaquette, then the contribution from four holes on 
this plaquette to the Gibbs weight of the partition func- 
tion will be just unity but the contribution of two paral- 
lel dimers arranged either vertically or horizontally will 



be 2z 



When 2z e u = 1, the holes become unstable 



to the formation of pairs of parallel nearest neighboring 
dimers and the result is a columnar dimer crystal. 

When C 2 = C 4 = (whereas C 6 remains positive), 
there is a mean-field tricritical point where the transition 
ultimately becomes continuous. Using the approximate 



Eqs. (|A43[) - (|A45[) . we can have an estimate for the tri- 
critical point. Solving the set of equations we arrive to 
the following estimate: 



z tr = 0.075 
u tr = 2.733 



(A49) 
(A50) 



Remarkably enough, these estimates are very close to 
the estimates from the grand-canonical simulations that 
are presented in Section |Vj However, we should be very 
cautious on taking these estimates too seriously since the 
terms in the expansion of the coefficients Co , C 2 , C4 , Cq 
in powers of z, generically have alternating signs with in- 
creasingly large constant coefficients which suggests that 
this expansion is not convergent. In any case, mean- 
field approximations, such as the one presented here, fail 
in two dimensions. The actual multicritical point has a 
more complex analytic structure, akin to a Kostcrlitz- 
Thouless transition, than suggested by these Landau- 
Ginzburg type arguments. On the other hand, the mul- 
ticritical point can be approached from the high dimcr 
density limit, where the transition is continuous and thus, 
an effective field theory description is possible. As we 
show in the following Sections IIV AlIVBl this multicriti- 
cal point is controlled by a marginally relevant operator 
and belongs to the Kostcrlitz-Thouless universality class. 
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APPENDIX B: DERIVATION OF THE SMA 
OSCILLATOR STRENGTH FUNCTIONS /(k). 

In this Appendix we present the details of the calcula- 
tions of the SMA oscillator strength functions /(k) dis- 
cussed in Section IVIII 



hole density operator commutes with the dimer density 
operator. Now, we will repeatedly use that 



,0- 



T2ct' : 



(B2) 



1. The fixed hole-density model 

a. Hole density excitations 



The only term of the Hamiltonian (|2.6j) which does not 
commute with the hole density operator is the hopping 
term for the holes. This term, in terms of destruction- 
creation Pauli operators <r^ , cr h can be written as: 



r t _ 



hole 



-f~hoie £ a h (r,)<r h+ (r fc )og + (r«)og (v jk ) (Bl) 

<ijk> 



At this point, we are not interested for the possible ori- 
entations of the dimers with respect to the holes, as the 



at the same position in real space. At any distance dif- 
ferent from zero, the commutator vanishes. For a given 
hole at a position R, there are four possible positions 
R' = R ± x ± y where it may move through the avail- 
able hopping term. By counting contributions from all 
these terms for every site of the lattice, we exactly take 
into account all the terms of the Hamiltonian including 
Hcrmitean conjugates. 

If R is the initial hole position and R + ro the final 
one, then ro has eight possible values: ro = ±i ± y. 
The first commutator can now be computed for any site 
R(The operator 7hoio(R, i"o) contains each of the above 
eight hopping terms). We have: 



[-t holc T holc (R, r ), a h (k)] = -t holc T holc (R, r ) (V lkR - e -' ik '( R + r »)) (B3) 

[o- h (-k), -t holc T holc (R 7 r )] = t h oieT hole (R, r ) (V kR - e lk '( R+r °)) (B4) 

[a /l (-k), [-t ho i e T hole (R,r ),cr ,l (k)]] = 2t holc 7i 10 i (R, r ) (1 - cos(k • r )) (B5) 
Thus, the oscillator strength /(k) can now be calculated: 

/(k) = ^([^(-k), [-t holc 7i lolc (R,ro),a' l (k)]])2t holc £(l-cos(k.r )) (B6) 

R.ro r 



If we set k = Qo+q where q is assumed to be small, then: 
For Q = (0,0), the above expression can be expanded: 

/(k) - 4t holc q 2 (B7) 
For Qo = (tt, 7r) we have: 

/(k) = 2t holc J2 [1 - (-ir-+^ cos(q • r )] (B8) 

ro 

fox + r 0y can take only the values and 2. Thus, 
(_iyo*+ra y — i f or a u ro ' s _ j n ^his wa y ; jf we expand 

around (7r, 7r), we have again: 

/(k) = 4t holc q 2 (B9) 

It is obvious that for any other k's, the oscillator strength 
/(k) cannot vanish. So, the upper bound to the excita- 
tion energy for the holes is: 

Ek — Eq < —77-r (BIO) 



b. Dimer density excitations 

In the same way as before, we may calculate the nu- 
merator for the case where we use the dimer density op- 
erator 0% instead of the hole density. Now, the operator 
does not commute with both the dimer-flip term and the 
hole-hopping terms. 

In the case of the dimer-flip term the commutator will 
give the following contribution^^: (the dimer-density op- 
erator is taken to be at a — x direction) 

/dimer-flip(k) = 8t £ [1 + COs(k • r )] (Bll) 
r=±y 

where ro = ±y for the case of a horizontal dimer. This 
term comes from the original quantum dimer model. As 
it was pointed out in Ref. l8rl it vanishes quadratically 
at Qo = (71", 7r) and Q\ — (0,7r). In particular, at Qo = 
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(-7T, 7r), with k = Qo + q, it vanishes as 

jdimer-flip(q) = 8tq 2 

In the case of hole-hopping terms which mix horizontal 



and vertical dimers r' Q — — Yq/2 = (±x ± y)/2 (where 
ro denotes the displacement vector for the hole in the 
(B12) considered move), we have: 



4(-k), [-i ho i e r« e (R, r ), at(k) 



= -t holc r h W(R,ro)f e - k -( R+ -/ 2 



ik-R 



2t holc T^l(R, r )(l - cos(k ■ r /2)) 



/hole(i)(k) = 2t hole ^(l - cos(k • r /2)) (B15) 

2. The fixed hole fugacity model 

Let's follow the same strategy as before (r = x,y): 



(B13) 
(B14) 



-t hoic T^ c (R, r ),a h (k) = t holc T^ c (R, r ) ( e 



— ■ik-R. _|_ g — 'ik-( 



( 2") 

where Tj^J denotes the resonance part of the Hamiltonian in Eq. (|2.8p . We also have, 



a h (-k), [-thoie^CR, r ), <7 ft (k)] ] = 2t holc T h ( o 2 1 ) c (R, r )(l + cos(k ■ r )) 
Finally, after adding the contribution of the hcrmitean conjugate part of the Hamiltonian, we have: 

/holo(2)(k) = 4t p 

airing 

(2 + cos(k a ) + cos(k y )) 
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